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Abstract 


During the last phase of the project, emphasis has changed to flexible space robotics, 
by mutual agreement between Dr. R. C. Montgomery, NASA Technical Officer, and the 
Principal Investigator. 

Significant advances have been achieved over the period covered by this report. Research 
has been concerned with two main subjects: 1) the maneuvering and control of freely floating 
flexible space robots and 2) the development of a theory for the motion of flexible multibody 
systems. Work on the first subject has resulted in two papers, both of them concerned with 
planar maneuvers. The first is concerned with the maneuvering and delivery of a payload 
to a certain point and in a certain orientation in space. The second is concerned with 
the docking maneuver with a target whose motion is not known a priori. Both papers will 
appear in the Journal of Guidance, Control, and Dynamics. The second subject is concerned 
with the development of hybrid (ordinary and partial) differential equations for the three- 
dimensional motion of flexible multibody systems, a subject of vital interest in flexible space 
robotics. The paper will appear in the Journal of Guidance, Control and Dynamics in an 
issue dedicated to the memory of Lawrence W. Taylor, Jr. 

Abstracts and copies of the papers are hereby included. 



1. Meirovitch, L. and Lim, S., “Maneuvering and Control of Flexible Space Robots,” NASA 
Workshop on Distributed Parameter Modeling and Control of Flexible Aerospace Systems , 
Williamsburg, VA, June 8-10, 1992. Also Journal of Guidance, Control , and Dynamics 
(in press). 

This paper is concerned with a flexible space robot capable of maneuvering payloads. 
The robot is assumed to consist of two hinge-connected flexible arms and a rigid end- 
effector holding a payload; the robot is mounted on a rigid platform floating in space. 
The equations of motion are nonlinear and of high order. Based on the assumption that 
the maneuvering motions are one order of magnitude larger than the elastic vibrations, a 
perturbation approach permits design of controls for the two types of motion separately. The 
rigid-body maneuvering is carried out open loop, but the elastic motions are controlled closed 
loop, by means of discrete-time linear quadratic regulator theory with prescribed degree of 
stability. A numerical example demonstrates the approach. In the example, the controls 
derived by the perturbation approach are applied to the original nonlinear system and errors 
are found to be relatively small. 

2. Chen, Y. and Meirovitch, L., “Control of a Flexible Space Robot Executing a Docking 
Maneuver,” AAS/AIAA Astrodynamics Conference , Victoria, B.C., Canada, August 16- 
19, 1993. Also Journal of Guidance , Control, and Dynamics (to appear). 

This paper is concerned with a flexible space robot executing a docking maneuver with a 
target whose motion is not known a priori. The dynamical equations of the space robot are 
first derived by means of Lagrange’s equations and then separated into two sets of equations 
suitable for rigid-body maneuver and vibration suppression control. For the rigid-body 
maneuver, on-line feedback tracking control is carried out by means of an algorithm based 
on Liapunov-like methodology and using on-line measurements of the target motion. For the 
vibration suppression, LQR feedback control in conjunction with disturbance compensation 
is carried out by means of piezoelectric sensor/actuator pairs dispersed along the flexible 
arms. Problems related to the digital implementation of the control algorithms, such as 
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the bursting phenomenon and system instability, are discussed and a modified discrete-time 
control scheme is developed. A numerical example demonstrates the control algorithms. 

3. Meirovitch, L. and Stemple, T. “Hybrid Equations of Motion for Flexible Multibody 
Systems Using Quasi-Coordinates,” AIAA Guidance, Navigation, and Control 
Conference , Monterey, CA, August 9-11, 1993. Also Journal of Guidance, Control, and 
dynamics - Issue dedicated to L. W. Taylor, Jr. (to appear). 

A variety of engineering systems, such as automobiles, aircraft, rotorcraft, robots, 
spacecraft, etc., can be modeled as flexible multibody systems. The individual flexible bodies 
are in general characterized by distributed parameters. In most earlier investigations they 
were approximated by some spatial discretization procedure, such as the classical Rayleigh- 
Ritz method or the finite element method. This paper presents a mathematical formulation 
for distributed-parameter multibody systems consisting of a set of hybrid (ordinary and 
partial) differential equations of motion in terms of quasi-coordinates. Moreover, the 
equations for the elastic motions include rotatory inertia and shear deformation effects. 
The hybrid set is cast in state form, thus making it suitable for control design. 
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MANEUVERING AND CONTROL OF FLEXIBLE SPACE ROBOTS* 


Leonard Meirovitch* and Seungchul Lim** 

Department of Engineering Science and Mechanics 
Virginia Polytechnic Institute and State University 
Blacksburg, VA 24061 

ABSTRACT 

This paper is concerned with a flexible space robot capable of maneuvering payloads. 
The robot is assumed to consist of two hinge- connected flexible arms and a rigid end-effector 
holding a payload; the robot is mounted on a rigid platform floating in space. The equations 
of motion are nonlinear and of high order. Based on the assumption that the maneuvering 
motions are one order of magnitude larger than the elastic vibrations, a perturbation ap- 
proach permits design of controls for the two types of motion separately. The rigid-body 
maneuvering is carried out open loop, but the elastic motions are controlled closed loop, by 
means of discrete-time linear quadratic regulator theory with prescribed degree of stability. 
A numerical example demonstrates the approach. In the example, the controls derived by 
the perturbation approach axe applied to the original nonlinear system and errors are found 
to be relatively small. 

1. INTRODUCTION 

A variety of space missions can be carried out effectively by space robots. These mis- 
sions include the collection of space debris, recovery of spacecraft stranded in a useless orbit, 
repair of malfunctioning spacecraft, construction of a space station in orbit and servicing the 
space station while in operation. To maximize the usefulness of the space robot, the manip- 
ulator ar ms should be reasonably long. On the other hand, because of weight limitations, 
they must be relatively light. To satisfy both requirements, the manipulator arms must be 

* Research supported by the NASA Research Grant NAG-1-225 monitored by Dr. R. C. 
Montgomery. The support is greatly appreciated. 

* University Distinguished Professor 
** Graduate Research Assistant 
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highly flexible. Hence, space robots share some of the dynamics and control technology with 
articulated space structures. 

Robotics has been an active research area for the past few decades, but applications 
have been concerned primarily with industrial robots, which axe ground based and tend to 
be very stiff and bulky. In contrast, space robots axe based on a floating platform and tend 
to be highly flexible. Hence, whereas industrial and space robots have a number of things 
in common, the differences axe significant. More recent investigations have been concerned 
with flexible industrial robots. 1-4 On the other hand, some investigations axe concerned 
with space robots consisting of rigid links. 5-7 Research on flexible 3pace robots has come 
to light only recently. 8,9 

This paper is concerned with a flexible space robot capable of maneuvering payloads. 
The robot is assumed to consist of two hinge-connected flexible arms and a rigid end-effector 
holding a payload; the robot is mounted on a rigid platform floating in space (Fig. 1). The 
platform is capable of translations and rotations, the flexible arms axe capable of rotations 
and elastic deformations and the end-effector/payload can undergo rotations relative to the 
connecting flexible arm. Based on a consistent kinematical synthesis, the motions of one body 
in the chain takes into consideration the motions of the preceding body in the chain. This 
permits the derivation of the equations of motion without the imposition of constraints. The 
equations of motion axe derived by the Lagrangian approach. The equations are nonlinear 
and of relatively high order. 

Ideally, the maneuvering of payloads should be carried out without exciting elastic vi- 
bration, which is not possible in general. However, the elastic motions tend to be small 
compared to the rigid-body maneuvering motions. Under such circumstances, a perturba- 
tion approach permits separation of the problem into a zero-order problem (in a perturbation 
theory sense) for the rigid-body maneuvering of the space robot and a first-order problem 
for the control of the elastic motions and the perturbations from the rigid-body motions. 
The maneuvering can be carried out open loop, but the elastic and rigid-body perturbations 
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axe controlled closed loop. 

The robot mission consists of carrying a payload over a prescribed trajectory and plac- 
ing it in a certain orientation relative to the inertial space. For planar motion, the end- 
effector/payload configuration is defined by three variables, two translations and one rota- 
tion. At the end of the mission, the vibration should be damped out, so that the robot can 
be regarded sis rigid at that time. Still, the rigid robot possesses six degrees of freedom, 
two translations of the platform and one rotation of each of the four bodies, including the 
platform. This implies that a kinematic redundancy exists. This redundancy can be used to 
optimize the robot trajectory 10 in the context of trajectory planning. A simpler approach is 
to remove the redundancy by imposing certain constraints on the robot trajectory, such as 
prescribing the motion of the platform . 11 Then, for a given end-eff ector/payload trajectory, 
the rigid-body maneuvering configuration of the robot can be obtained by means of inverse 
kinematics. Finally, the forces and torques required for the robot trajectory realization are 
obtained from the zero-order equations by means of inverse dynamics. 

The first-order equations for the elastic motions and the perturbations in the rigid-body 
maneuvering motions axe linear, but of high order, time-varying and they are subjected 
to persistent disturbances. The persistent disturbances arise from the zero-order solution, 
and hence axe known; they axe treated by means of feedforward control. All other distur- 
bances are controlled closed loop, with the feedback control being designed by means of 
discrete- time linear quadratic regulator (LQR) theory with prescribed degree of stability. A 
numerical example demonstrates the approach. In the example, the controls derived by the 
perturbation approach are applied to the original nonlinear system and the errors in the end 
effector/payload configuration were found to be relatively small during the maneuver and to 
vanish soon after the termination of the maneuver. 

2. A CONSISTENT KINEMATICAL SYNTHESIS 

To describe the motion of the space robot, it is convenient to adopt a consistent kine- 
matical synthesis whereby the system is regarded as a chain of articulated flexible bodies and 
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the motion of one body is defined with due consideration to the motion of the preceeding 
body in the chain. Figure 1 shows the mathematical model of a planar space robot. The 
robot consists of a rigid platform (Body 1), two hinge-connected flexible arms (Bodies 2 and 
3) and a rigid end-effector holding the payload (Body 4). The various motions are referred 
to a set of inertial axes and sets of body axes to be defined shortly. 

The object is to derive the system equations of motion, which can be done by means of 
Lagrange’s equations in terms of quasi-coordinates. 12 Because in the case at hand the motion 
is planar, it is more expedient to use the standard Lagrange’s equations. This requires the 
kinetic energy, potential energy and virtual work. The kinetic energy, in turn, requires the 
velocity of a typical point in each of the bodies. 

The position of a nominal point on the platform is given by 

Ri = Ro + ri (1) 

rrs 

where Ro = [X Y] is the position vector of the origin 0\ of the body axes xi,yi (Fig. 1) 
relative to the inertial axes X, Y and in terms of A, Y components and ri = [xi y{\ is the 
position vector of the nominal point on the platform relative to the body axes x\ , y\ and in 
terms of xi, yi components. The velocity vector of a point on the platform can be expressed 
in terms of x\,y\ components as follows: 

V 1 = C 1 Ro + u l r 1 (2) 


where 


(3) 

is a matrix of direction cosines between axes xi,yi, and X, Y, in which a6\ = sin#i, cB\ = 
cos 6 \ , 





(4) 


is the velocity vector of 0\ in terms of A, Y components and 


*1 


0 -0i 

(9i 0 


(5) 
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The second body is flexible, so that we must resolve the question of body axes. We 
define the body axes 22, yi as a set of axes with the origin at the hinge O2 and embedded in 
the undeformed body such that 12 is tangent to the body at O 2 (Fig. 2). Then, we define 
the motion of axes 22, y 2 as the rigid-body motion of Body 2 and measure the elastic motion 
relative to 22, y2- Hence, the velocity of a point on Body 2 (first flexible arm) in terms of 
S 2 ,y2 components is 

V2 =£2— lVi ( O2 ) + &2 (r 2 + U2) + u 2re i 

=C2Rq + C2-iwiri ( O2 ) + (r2 + U2) 4- u 2re } (6) 

where C2-1 and C2 are matrices similar to C\, Eq. (3), except that 6 \ is replaned by $2— 0 \ and 
82, respectively, £>2 has the same structure as £>1 but with 62 replacing 9 \, ri (0 2) = [di A-i] , 

rp rp 

r2 = [X2 0 ] , U2 = [0 U2] and u 2re i = [0 U2], in which U2 = U2 (12 ,t) and U2 = U2 (22 ,t) 
are the elastic displacement and velocity, respectively. 

Using the analogy with Body 2, the velocity of a point on Body 3 (second flexible arm) 
in terms of 23, y3 components can be shown to be 

V3 —Cz- 2^2 (-£2) +£3 (r3 + U3) + Uj r gj 

=C3Ro + Cs-i^iri (O2) + £3-2 \&2 [ r 2 (^2) + u 2 (-^2) 0) ^2rel C^2» 0} 

+ W3(r3 + u 3 ) + u 3rel (7) 

The fourth body consists of the end-effector and payload combined, and is treated as 
rigid. Following the established pattern, the velocity of a point on Body 4 in terms of 24, y4 
components is 

V 4 =C f 4_3V3 (£3) +£>41*4 

=(^4^ + ( O2 ) + ^4-2 \&2 [ r 2 (-^2) + u 2 (-^2»0] + ^2rel 

+ C 4 — 3 {^3 [r3 (-^ 3 ) + U 3 (£ 3 , 0 ] + ^3rel (^3, *)} + ^4^4 (8) 

The consistent kinematical synthesis just described permits the formulation of the equa- 
tions of motion for the whole system without invoking constraint equations. Such constraint 
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equations must be used to eliminate redundant coordinates in a formulation in which equa- 
tions of motion are derived separately for each body. 

3. SPATIAL DISCRETIZATION OF THE FLEXIBLE ARMS 

The velocity expressions derived in Sec. 2 involved rigid-body motions depending on 
time alone and elastic motions depending on the spatial position and time. Equations of 
motion based on such formulations are hybrid, in the sense that the equations for the rigid- 
body motions are ordinary differential equations and the ones for the elastic motions are 
partial differential equations. Designing maneuvers and controls on the basis of hybrid 
differential equations is likely to cause serious difficulties, so that the only viable alternative 
is to transform the hybrid system into one consisting of ordinary differential equations alone. 
This amounts to discretization in space of the elastic displacements, which can be done by 
means of series expansions. Assuming that the flexible arms act as beams in bending, the 
elastic displacements can be expanded in the series 

u i (*»>*) = i ( x ») Vij (<) = 4>l » * = 2,3 (9) 

i-i 

where <f>(j (x*) axe admissible functions , often referred to as shape functions , and rjij ( t ) 
axe generalized coordinates (i = 2,3; j = 1,2, ...,rii); <j>{ and are corresponding rij- 
dimensional vectors. 

The question arises as to the nature of the admissible functions. Clearly, the object is to 
approximate the displacements with as few terms in the series as possible. This is not a new 
problem in structural dynamics, and the very same subject has been investigated recently in 
Ref. 13, in which a new class of functions, referred to as quasi-comparison functions^ has been 
introduced. Quasi- comparison functions are defined as linear combinations of admissible 
functions capable of satisfying the boundary conditions of the elastic member. As shown in 
Fig. 2, the beam is tangent to axis x, at 0{ (t = 2,3). Hence, the admissible functions must 
be zero and their slope must be zero at Xj = 0. At x, = the displacement, slope, bending 
moment and shearing force axe generally nonzero. Quasi-comparison functions are linear 
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combinations of functions possessing these characteristics. Admissible functions from a single 
family of functions do not possess the characteristics, but admissible functions from several 
suitable families can be combined to obtain them. In the case at hand, quasi-comparison 
functions can be obtained in the form of suitable linear combinations of clamped-free and 
clamped-clamped shape functions. 

4. LAGRANGE’S EQUATIONS 

Before we can derive Lagrange’s equations, we must produce expressions for the ki- 
netic energy, potential energy and virtual work. To this end, and following the spa- 
tial discretization indicated by Eqs. (9), we introduce the configuration vector q(t) = 

p 

[X(t) Y(t) 6i(t) 0 2 (t) 0 3 (<) 04 (<) rffit) Vz{t)} so that the velocity vectors, Eqs. (2), (6)- 
(8), can be written in the compact form 

V i = Diq, i— 1,2, 3, 4 (10) 


where 


Di = 

D 2 = 


cB 1 

s6\ 

-yi 0 ... o r l 



. —s9i 

cB\ 

XI 0 ... 0 r . 



c6 2 

30 2 

di3 ( 0 2 — 9\) — hic(0 2 — 9i ) 

0 

0 

0 

Cr 

1 

—30 2 

cB 2 

dic(0 2 — 9i) + his ( 0 2 — #i) 

x 2 0 0 <t> 


0 T ' 

nT 


( 11 ) 


Then, the kinetic energy is simply 

T = \ £ / VfVrfmi = l -q T Mq ( 12 ) 

t t=1 Jrni * 

where 



is the mass matrix. Typical entries in the mass matrix are 


mn = m, mi 2 = 0, m 13 = - (m 2 + m 3 + m 4 ) {hic&i + dis^i) 
m 14 = - ^2 + + >™ 4 ) <f>2 (^2)j V2 c &2 - [^2 + ( m 3 + m 4) L 2 ] S0 2 
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(14) 


mis = - 


-rT 


<t> 3 + rrn<j ) 3 (L 3 ) 


3 & 3 


m.22 = m, m23 = — ( m 2 + m3 + T774) (h^sOi — d\cB\) 


m28 = + rrn4$ (£ 3 ) c#3 


m 8 8 = L , <t> 3 <j) 3 dm 3 + m 4 <f> 3 (L 3 )<j> 3 (L 3 ) 

J Body 3 


in which 


m = ^ m,, = / 4>idmi, i = 2,3, S'j = / xidrrii , t = 1,2,3, 4 

The potential energy, assumed to be entirely due to bending, has the form 

V = Zj Q 2 Eh u" (x 2 , t) 2 <£32 + ^ jf EI 3 [«3 (*3, *)] 2 dx 3 = ^q T Kq 


(15) 


(16) 


in which £/, (1 = 2,3) are bending stiffnesses and primes denote spatial derivatives. More- 
over, 

K = block-diag [0 K 2 K 3 ] (17) 


is the stiffness matrix, where 

Ki = / L< Elitftfdxi, i = 2,3 (18) 

axe stiffness matrices for the flexible arms. 

Next, we propose to derive the virtual work expression. To this end, we must specify first 
the actuators to be used. There are three actuators acting on the platform, two thrusters 
F z \ and F y i acting at 0\ in directions aligned with the body axes and a torquer M\. Three 
other torquers M2, M3 and M4 axe located at the hinges 02, 03 and O4, respectively, the first 
acting on the platform and first arm, the second acting on the first and second arm and the 
third acting on the second arm and end-effector. Moreover, there are torquers M$, M$ y M7 
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and Mg acting at xi = X 2 / 3 , X 2 = 2 X 2 / 3 , X 3 = X 3/3 and X 4 = 2 X 3 / 3 , respectively. In view 
of this, the virtual work can be written as follows: 


SW =F t 1 (cos 9\8X 4- sin 9\ 8 Y) + F v \ (— sinfli&Y + cos0i<5Y) + M\89\ 
+ M 28 ( 02 — #i) + Mz 8^2 + M± 8 \j)± 4- M $8 \$2 + 4>2 (£ 2 / 3 ) t7 2 
+ M 6 8 [0 2 + (2X2/3) rj 2 ] + M 7 8 [0 3 + (X 3 /3) 173 ] 

+ Mg [^3 4- 4>3 ( 2 X 3 / 3 ) T 73 


(19) 


where 8 X, 8 Y , . . . axe virtual displacements. Moreover, denoting the angles between the two 
arms and between the second arm and the end-effector by 


03 =03 — 02 ~ 


04 =04 ~ 03 ~ 


du2 

3x2 

duz 

8 x 3 


12=1.2 


*3 = ^3 


— 03 — 02 ~ 4 * 2 ' (^ 2 ) V 2 

= 6, -Bs -<t>f (Ls)sis 


( 20 ) 


we can write 


8+3 = 89 3 - 80 2 - 4>2 (L 2 ) 8v 2 , *04 - 80 3 - <f>? (X 3 ) ^3 ( 21 ) 

Inserting Eqs. (21) into Eq. (19), we can express the virtual work in terms of generalized 
forces and generalized virtual displacements in the form 

SW = Q r £q (22) 

2" 1 

where Q = [f x F y ©i 0 2 ©3 ©4 Nf] is the generalized force vector, in which 

F x = F t 1 cos 9\ — F y i sin 9\, Fy = F x 1 sin 9\ + F y \ cos 6\ 

0i = Mi — M2, @2 = M2 — M3 + Mj + M5 + Ms 

03 = M3 — M4 + Mj 4- Mg, 04 = M4 (23) 

N 2 = -M 3 d>2 (X 2 ) + M 5 d>2 (X2/3) + Ms 02 (2X2/3) 

Na = -M 4 0<, (X 3 ) + M 7 0' 3 (X 3 /3) + m 8 0' 3 ( 2X3/3) 
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X* 

and <5q = ^8X 8Y 89\ 89 2 69$ ^4 $12 is the generalized virtual displacement vector. 

Equations (23) express the generalized forces and torques in terms of the actual actuator 
forces and torques and can be written in the compact form 

Q = EF (24) 

rr\ 

where F = [F x i F y \ M\ M 2 . . . M%\ is the actual control vector and 
E = E{6i) = 

ci —si 000 0 0 0 0 0 

si ci 000 0 0 0 0 0 

001-10 0 0 0 0 0 

0 0 0 1 -1 0 1 1 0 0 
00001 -1 0 0 1 1 

00000 1 0 0 0 0 

0000 -<A' 2 (L 2 ) 0 4 >\ ( y) 0 0 

0000 0 -< i >' 3 ( L 3 ) 0 0 ^( y ) ^( x ) 

(25) 


where si = sin^i, ci = cos#i. Note that E is a time- varying coefficient matrix, because 9\ 


varies with time. 

Lagrange’s equations can be expressed in the general symbolic vector form 

£(dT\ = ar ev = 

dt \dq) dq 8q w ^ 1 

Observing that M = M (q), we can write 

3T d (dT\ 1# . • 

a* =A/q> q+ q 

dT 1 , t 9M. dV 

K = 2‘ t K q ’-^ = Kq (2?) 

Inserting Eqs. (27) into Eq. (26), we obtain Lagrange’s equations in the more explicit form 

Mq+(Af-iq r ^)q + tfq=Q (28) 
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in which 


• 6 +2*<9M. 

M = 2_s <7i> 

U d( n 


* dM_ 

dq 


q^dM/dqi 

q T dM/dq2 

m q T dM/dq6+2n. 


5. A PERTURBATION APPROACH TO THE CONTROL DESIGN 


(29) 


Equation (29) represents a high-order system of nonlinear differential equations, and is 
not very suitable for control design. Hence, an approach capable of coping with the problems 
of high- dimensionality and nonlinearity is highly desirable. Such an approach must be based 
on the physics of the problem. The ideal maneuver is that in which the robot acts as if 
its arms were rigid. In reality, the arms are flexible, so that some elastic vibration is likely 
to take plane. It is reasonable to assume, however, that the elastic motions are one order 
of magnitude smaller that the maneuvering motions. This permits treatment of the elastic 
motions as perturbations on the maneuvering motions. In turn, the elastic perturbations give 
rise to perturbations in the “rigid-body” maneuvering motions. This suggests a perturbation 
approach, whereby the problem is separated into a zero-order problem for the “rigid- body” 
maneuvering of the payload and a first-order problem for the control of the elastic motions 
and the perturbations in the rigid-body maneuvering motions. The zero-order problem is 
nonlinear, albeit of relatively low dimension. It can be solved independently and the control 
can be open loop. On the other hand, the first-order problem is linear, but of relatively high 
dimension. It is affected by the solution to the zero-order problem, where the effect is in the 
form of time-varying coefficients and persistent disturbances. The control for the first-order 
problem is to be closed loop. 

We consider a first-order perturbation solution characterized by 


q = q 0 +qi, Q = Qo + Ql 


(30) 


where the subscripts 0 and 1 denote zero-order and first-order quantities, with the zero-order 
quantities being one order of magnitude larger than the first-order ones. Inserting Eqs. (30) 
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into Eq. (28), separating quantities of different orders of magnitude and ignoring terms of 
order two and higher, we obtain the equation for the zero-order problem 


Moqo + ^ M v — qo — Qo — FJqFq 


where qo = [Jfo Vo #10 #20 #30 #40 0 r 0 r ] , Qo = [Fxo Fyo ©10 ©20 ©30 ©40 0 T 0 r ] are 
zero-order displacement and generalized control vectors, Eq = E (# 10 ) is the matrix E, 
Eq. (25), evaluated at 9\ = <?io, Fq = [F x 0 F v q M\q M 20 ... M8o] r and 


„ w/ \ „ \ 9 M. dM. dM ^ 1 

Mq = M (qo) , M„ = - — qo -5 — qo • • • t, qo 

. Oq\ Oq 2 O l ?6+2n J q=qo 

Moreover, we obtain the equation for the first-order problem 


(32a, b) 


M 0 qi + (M v + M'~ M?) q! 4- yM a + M vv - -M' vv + K J qi = Qi + Q d (33) 

where qi = [Xi Yi B\\ #21 #31 #41 I 2 V 3 , Qi = [-Fjn Fy\ ©n ©21 ©31 ©41 N^j 
are first-order displacement and generalized control vectors, Qj =|o 000000 Fj 2 F^l 
is a persistent disturbance vector and 


, r dM . dM _ dM _ 1 

M a = . qo 0 qo • • • * qo 

Oq\ Oq 2 o?6+2n J q=qo 

(34a) 

6+2,1 dM 

* = £ 57 <*>; 

0( Jj q=qo 

(346) 

6+2n64-2n ^2^ 

M„t,qi = V) ]T) ~ ~ ?iJt?Ojqo 

y_l jt = i dqjdqk q=q 0 

(34c) 

>r < 6+2,1 a 2 M 

^««qi = qo L x n x n ?i*qo 

jfc=l OQvqk q=q 0 

(34d) 


From Eqs. (24) and (25), however, we can write 

Qi = EoF 1 + E1F0 = F0F1 + FJqi 

where 

p dE 0 

off 1 ex=«io 
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Moreover, the matrix Fq has the entries 

foil = - (F z 10 sin 0io + F y 10 cos 0io) 

F 021 = Fzio cos ^10 - F y 10 sin 0io (37) 

Foij = 0, i = 3,4, . . . , y + n 2 + n 3 \ j = 2,3,... ,6 + n 2 + n 3 
In view of this, the equation for the first-order problem, can be rewritten as 

M 0 qi + ( M v + M' - A/J) qi + {Mi + M vv - 2^00 + K ~ F^j qi = EqFi + (38) 

6. TRAJECTORY PLANNING 

The mission consists of delivering the payload to a certain point in space and placing it in 
a certain orientation. For planar motion, the final payload configuration is defined by three 
variables, two translations and one rotation. The trajectory planning, designed to realize this 
final configuration, will be carried out as if the robot system were rigid, with the expectation 
that all elastic motions and perturbations in the rigid-body maneuvering motions will be 
annihilated by the end of the maneuver. The rigid-body motion of the robot is described by 
the zero-order problem and it consists of six components, two translations of the platform and 
one rotation of each of the four bodies. This implies that a kinematical redundancy exists, as 
there is an infinity of ways a six-dimensional configuration can generate a three-dimensional 
trajectory. This redundancy can be removed by controlling surplus variables, perhaps in an 
optimal fashion. In this study, we prescribe three of the configuration variables, such as the 
translations and rotation of the platform. Under these circumstances, the rigid space robot 
can be treated as a nonredundant manipulator. 

Next, we denote the end-effector configuration by X £, so that from kinematics we can 
write 

X £ = f(q 0 ) (39) 

where f is a three-dimensional vector function. Differentiating Eq. (39) with respect to time, 
we obtain 

X E = J{ qo)qo (40) 
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where 


J (qo) = [df/d q 0 ] (41) 

is the 3x6 Jacobian matrix. Introducing the notation 

r ( -j T 

qo = q 5 ! qlr (42) 

rp rp 

where = [Xo Yq #io] and q m = [#20 #30 # 4 o] are the controlled platform configura- 
tion vector and the open-loop controlled manipulator configuration vector, respectively, and 
partitioning the Jacobian matrix accordingly, or 

J = Jx • Jm (43) 

Eq. (40) can be rewritten as 

X E = Jsqs + Jmqm (44) 

Then, on the assumption that q£ is prescribed and for a given end-effector trajectory X E , 
we can determine the manipulator velocity vector from 

qj* = Jm {Xe - Js*s) (45) 

The end-effector trajectory was taken in the form of a sinusoidal function so as to prevent 
excessive vibration. Finally, with qo given, we can obtain the required open- loop control F o 
by inverse dynamics, which amounts to using Eq. (31). 

7. FEEDBACK CONTROL OF THE ELASTIC MOTIONS AND RIGID-BODY 
PERTURBATIONS 

The elastic motions and the perturbations in the rigid-body maneuvering motions are 
governed by the equation defining the first-order problem, Eq. (38). The persistent distur- 
bances are controlled open loop and all other disturbances are controlled closed loop. To 
this end, we express the control vector in the form 

Fi = Fi„ + F lc (46) 
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where the subscripts o and c indicate open loop and closed loop, respectively. Recognizing 
that Eq is a rectangular matrix, the open-loop control can be written as 



Flo — FoQd 

(47) 

in which 




4 = (E^E o y' E% 

(48) 

is the psuedo-in verse of Eq. 




For the closed-loop control, we consider LQR control, which requires recasting the equa- 
tions of motion in state form. Adjoining the identity qi = qj, the state equations can be 
expressed as 

x(i) = A(i)x(i) + B(t)E 0 u c (t) + B(t)Dd(t) (49) 

X* 

where x = is the state vector, u c = Fi c is the control vector, d = is the 

disturbance vector and 

, r o i 

A - -M 0 -‘ (M„ + M„„ - | Ml, + K- F$) -M 0 -> [M. + M< - A/„ r 

B= [ Z> = (/-£o4) (506, c) 

axe coefficient matrices. It should be noted here that, if the matrix Eq is not square, the 
matrix D is not zero, so that the open-loop control does not annihilate the persistent distur- 
bances completely. As the number of actuators approaches the number of degrees of freedom 
of the system, the matrix Eq tends to become square. When the number of actuators co- 
incides with the number of degrees of freedom the matrix Eq is square, in which case the 
pseudo-inverse becomes an exact inverse and the matrix D reduces to zero. 

The state equations, Eq. (49), possess time- varying coefficients and axe subject to residual 
persistent disturbances. Due to difficulties in treating such systems in continuous time, we 
propose to discretize the state equations in time. Following the usual steps, 14 the state 
equations in discrete time can be shown to be 

Xfc+i = $ifcXjfe + T k E Q k\i c k + TjtDjtdjfc, A: = 0,1,... (51) 
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where 


x* =x(fcT), u ck = u c (kT), d t =d (JcT) y k = 0.1,... 

=exp A k T, T k = (exp A*r - I) A* 1 B kf k = 0, 1, . . . 
£ 0 * =E 0 (kT), D k = D (kT) , fc = 0,1,... 


(52) 


in which T is the sampling period and 

A* = A(fcT), B k = B(fcT) (53) 

In view of the above discussion, we assume that the effect of the persistent disturbances has 
been reduced drastically by the feedforward control, and design the feedback control in its 
absence. This design is according to a discrete-time LQR with prescribed degree of stability. 
To this end, we consider the performance measure 

J = x^P N x N + 53 e2ak ( x fc QjfeXfc + Udt-Rjfcii c *) (54) 

k = o 

where Pfj and Q k are symmetric positive semidefinite matrices, R k is a symmetric positive 
definite matrix, cc is a nonnegative constant defining the degree of stability and NT is the 
final sampling time. 

The optimization process using the performance measure given by Eq. (54) can be re- 
duced to a standard discrete-time LQR. form by means of the transformation 

x k = e ak x ki u ck = e“*u c jb, Pn = e~ 2aN P n (55a, b, c) 

Multiplying Eqs. (51) through by using Eqs. (55a, b) and ignoring the small perturbing 

term, we obtain the new state equations 

xjk+i = e“ (* k x k + 'y k Eo k u ck ) , k = 0, 1, . . . , N - 1 (56) 

Similarly, inserting Eqs. (55) into Eq. (54), we obtain the new performance measure 

J = XffpN^N + 13 {x-kQkX-k + u^ k R k u ck j (57) 

Jfc=0 
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It can be shown that the optimal control law has the form 14 


u c * = G k x k , k = 0, 1, . . . , N - 1 (58) 

where G k are gain matrices obtained from the discrete-time Riccati equations 

Gff-i = ~{e 2a EQ^^iT^f^ipN+l-iTN-iEo^-i+RN-i) 

i = 1 , 2 , . . . , N\ P = e~ 2 aN P N (59 a) 

rjy 

Ptf-i =e 2a ($N-i + ^N-i^ 0 ,N-iGN-i) P N+l-i ($N-i + ? N-iEo t fi[-iG fif-i] 

+ G%_iR N -iG N -i + Q N - iy i = 1 , 2 , . . . , AT; P N = e~ 2 aN P N (596) 

A A 

Equations (59a) and (59b) are evaluated alternately for Pn-u Gn- 2 , Pn-2 , • • • i Go, 

given the final value of P//. 

Inserting the control law, Eqs. (58), into Eqs. (56), we obtain the closed-loop transformed 
state equations 

x fc+ i = e“ (* fc + r*PoiG fc ) x fc , k = 0, 1 , . . . (60) 

Then, recalling Eq. (55a) and restoring the persistent disturbance term, the closed-loop state 
equations for the original system can be written in the form 

x* + i = ($* + T k E ok G k )x k + I^d*, A: = 0,1,... (61) 

8 . NUMERICAL EXAMPLE 

The example involves the flexible space robot shown in Fig. 1. Numerical values for the 
system parameters are as follows: 

L\ — 1 m, d\ — 0.5 m, I>2 = L3 = 5m., L4 = 1.66m 

mi = 10 kg, m 2 = m 3 = 1kg, m 4 = 0.1 kg 

J\ = 20 kgm 2 , J 2 = 3 kgm 2 , EI 2 = EI 3 = 122.28 Nm 2 
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The quasi-comparison functions for the flexible arm were chosen as a linear combination of 
clamped-free and clamped-clamped shape functions. Both families of shape functions have 
the functional form 

<f>i = - 4 = [cosh Xix/L — cos \{x/ L — ai (sinh A,x/L — sin Xix/L )] , t = 1,2, . . . ,n 
yjL 

The values of A, and cr, for each f ami ly are given in Table 1. They correspond to two 
clamped-free and three clamped-clamped shape functions, for a total of n = 5 for each 
flexible arm. 

The initial and final end-effector positions axe defined by 

Xi = 9.757 m, Y; = 1.914 m, 0 Ai = 0 rad 
Xf = 5.000 m, Yf = 1.914 m, 6±f = — tt/ 2 rad 

and we note that the path from the initial to the final position represents a straight-line 
translation, while the orientation undergoes a 90° change. In terms of time, the translational 
and rotational accelerations represent one-cycle sinusoidal curves. 

The maneuver time is tf — 2.5 s. The zero-order actuator forces and torques to carry 
out the maneuver are shown in Fig. 3. 

The control of the elastic motions and the perturbations in the rigid-body motions was 
extended to t = 4 s. Note that for 2.5 s < t < 4 s the system is time- invariant, during 
which time the control gains can be regarded as constant. The weighting matrices in the 
performance measure axe 

Qk = 10/, R k = /, Pn = 10/ 

The degree of stability constant is a — 0.1. Moreover, the samping period is T = 0.01 s and 
the number of time increments is N = 350. 

Time-lapse pictures of the uncontrolled and controlled robot configuration axe shown 
in Figs. 4a and 4b, respectively, at the instants 0, 1, 1.5 and 2.5 s. Figures 5 and 6 show 
time histories of the errors in the end-effector position. The discrete-time open-loop and 
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closed-loop poles for a = 0.01 axe given in Tables 2 and 3. For comparison, Fig. 7 shows the 
time history of the errors and Table 4 gives the closed-loop poles for a = 1. 

It should be pointed out that the actuator dynamics was also included in the computer 
simulation, but the effect turned out to be small. 11 

9. CONCLUSIONS 

An orderly kinematic synthesis in conjunction with the Lagrangian approach permits 
the derivation of the equations of motion for an articulated multibody system, such as 
those describing the dynamical behavior of a flexible space robot, without the imposition 
of constraints. The equations are nonlinear and of relatively high order. A perturbation 
approach permits the separation of the problem into a zero-order problem (in a perturbation 
sense) for the rigid-body maneuvering of the space robot and a first-order problem for the 
control of the elastic motions and the perturbations from the rigid-body motions. The 
robot mission consists of carrying a payload over a prescribed trajectory and placing it in 
a certain orientation relative to the inertial space. This represents the zero-order problem 
and the control can be carried out open loop. The first-order equations defining the first- 
order problem (in a perturbation sense) are linear, time-varying, of high-order and subject 
to persistent disturbances. The persistent disturbances are treated by means of feedforward 
control. All other disturbances axe controlled closed loop, with the feedback control being 
designed by means of discrete-time LQR theory with prescribed degree of stability. In a 
numerical example, the controls derived by the perturbation approach are found to work 
satisfactorily when applied to the original nonlinear system. 
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Table 1. Shape Function Coefficients 


i 

A, 

<?• 

i 

1.8751 

0.7341 

2 

4.6941 

1.0185 

3 

7.8548 

0.9992 

4 

10.9955 

1.0000 

5 

14.1372 

1.0000 


Table 2. Discrete-Time Open-Loop Poles 


No. 

Pole Location 

Mag. 

No. 

Pole Location 

Mag. 

1,2 

-0.840 ± 0.543i 

1.000 

17,18 

0.991 ± 0.135i 

1.000 

3,4 

-0.778 ± 0.629i 

1.000 

19,20 

0.994 ± 0.107i 

1.000 

5,6 

-0.700 ± 0.714i 

1.000 

21,22 

1.000 

1.000 

7,8 

-0.690 ± 0.724i 

1.000 

23,24 

1.000 

1.000 

9,10 

0.586 ± 0.810i 

1.000 

25,26 

1.000 

1.000 

11,12 

0.629 ± 0.778i 

1.000 

27,28 

1.000 

1.000 

13,14 

0.902 ± 0.431i 

1.000 

29,30 

1.000 

1.000 

15,16 

0.921 ± 0.390i 

1.000 

31,32 

1.000 

1.000 
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Table 3. Discrete-Time Closed- Loop Poles for a = 0.1 


No. 

Pole Location 



Pole Location 

Mag. 

1,2 

-0.169 ± 0.546i 

0.572 

18,19 

0.803 ± 0.976 x lO^i 

0.809 

3 

0.493 x 10" 2 

0.005 

20 

0.805 

0.805 

4 


0.012 

21 

0.807 

0.807 

5 

0.125 



0.814 ±0.362 x 10 -2 i 

0.814 

6 

0.204 



0.817 

0.817 

7,8 

0.302 ±0.148i 

0.336 

26 

0.817 

0.817 

9,10 

0.454 ± 0.493i 

0.670 

27 

0.819 

0.819 

EB9 

0.468 ± 0.323i 

0.569 

28,29 

0.821 ± 0.366 x 10 -2 i 

0.821 

12,13 

0.536 ±0.500i 

0.733 

30 

0.822 

0.822 

15,16 

0.749 ±0.860 x 10 _1 i 



0.822 

0.822 

17 

0.792 

0.792 

32 

0.827 

0.827 


Table 4. Discrete-Time Closed- Loop Poles for a = 1 


No. 

Pole Location 

Mag. 

No. 

Pole Location 

Mag. 

1 

-0.566 



0.139 ± 0.844 x 10" 2 i 

0.139 

2,3 

-0.160 ±0.186i 

0.246 

19,20 

0.150 ± 0.022i 

0.152 

4,5 

-0.109 ± 0.275i 

0.296 

21,22 

0.187 ± 0 . 1 45i 

0.236 

6,7 

0.062 ± 0.088i 

0.108 

23,24 

0.198 ± 0.288 x 10 _1 i 

0.200 

8 

-0.177 x lO" 1 

0.018 

25 

0.251 

0.251 

9,10 

0.779 x 10" 2 ± 0.209t 

0.209 

26,27 

0.252 ± 0.180i 

0.310 

11,12 

0.072 ±0.088i 

0.114 

28,29 

0.279 ± 0.490i 

0.564 

13,14 

0.118 ±0.016i 

0.119 

30,31 

0.328 ± 0.148i 

0.360 

15,16 

' 0.132 ±0.920 x 10 -2 i 

0.132 

32 

0.430 

0.430 
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CONTROL OF A FLEXIBLE SPACE ROBOT 
EXECUTING A DOCKING MANEUVER* 


Y. Chen* and L. Meirovitch** 

Department of Engineering Science and Mechanics 
Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 24061 

ABSTRACT 

This paper is concerned with a flexible space robot executing a docking maneuver with a 
target whose motion is not known a priori. The dynamical equations of the space robot axe 
first derived by means of Lagrange’s equations and then separated into two sets of equations 
suitable for rigid-body maneuver and vibration suppression control. For the rigid-body 
maneuver, on-line feedback tracking control is carried out by means of an algorithm based 
on Liapunov-like methodology and using on-line measurements of the target motion. For the 
vibration suppression, LQR feedback control in conjunction with disturbance compensation 
is carried out by means of piezoelectric sensor/actuator pairs dispersed along the flexible 
arms. Problems related to the digital implementation of the control algorithms, such as 
the bursting phenomenon and system instability, are discussed and a modified discrete-time 
control scheme is developed. A numerical example demonstrates the control algorithms. 

1. INTRODUCTION 

One of the functions of a space robot is to deliver payloads accurately and smoothly to a 
moving target. An example of such a space robot is shown in Fig. 1. The robot consists of a 
rigid base, two flexible arms attached to the base in series and an end-effector/payload. To 
carry out the mission described, the space robot must have its own control system enabling 
the platform to translate and rotate and its arms to rotate. In this paper, the target motion 
is assumed not to be known a priori, so that the control permitting the space robot to execute 

t Supported by the AFOSR Research Grant F49620-89-C-0045 monitored by Spencer 
T. Wu and by the NASA Research Grant NAG-1-225 monitored by Raymond C. 
Montgomery. 

* Graduate Research Assistant. 

** University Distinguished Professor, Fellow AIAA. 
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the docking maneuver must be based on on-line measurements. 

The equations governing the behavior of space robots are nonlinear and can be expressed 
in the general form of the state equation 

x = f(x,u) (la) 

and the output equation 

y = g(x) (16) 

where x is the state vector, u is the control force vector and y is the output vector, usually 
defined as the position and orientation variables of the end-effector. The target output vector 
y t is defined as the position and orientation variables of the target. We can then define the 
error vector as 

e = yt - y (2) 

The problem reduces to that of designing a control law u(t) so that e and its time derivative 
e are driven to zero. 

There are two significant differences between industrial robots in current use and space 
robots considered here. In the first place, industrial robots are mounted on a fixed base, 
whereas space robots axe mounted on space platforms capable of translations and rotations. 
The second significant difference is that space robots must be very light, and hence very 
flexible, unlike industrial robots characterized by very bulky and stiff arms. The flexibility 
of the robot arms causes elastic vibration, which tends to affect adversely the performance of 
the end-effector. Both a floating platform and flexibility axe being considered in this paper. 

In the case of space-based robots, research has been carried out on the assumption that 
the platform floats freely, 1-6 i.e., that there are no external forces and torques acting on the 
system, which implies that the system lineax and angular momentum axe conserved. For 
a space robot tracking a moving target, it is unrealistic to make such an assumption, so 
that algorithms concerned with free-floating space robots are not applicable to the problem 
considered here. 
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The most commonly used approach to robotics can be described as follows: first, inverse 
kinematics is performed to obtain the desired robot configuration trajectory q^(t) from the 
desired end-effector trajectory y <j(t). Then, using the system equations of motion, inverse 
dynamics is performed to obtain the control force realizing ^ the target motion is 

known a priori, the end-effector’s trajectory, as well as the robot trajectory, can be deter- 
mined by an off-line planning algorithm. For a kinematically- redundant robot, such as the 
one considered here, the robot redundancy can be used to achieve optimality. 7 

If the target motion is not known a priori, planning is impossible. Even when the target 
motion is known, it is very likely that some unexpected disturbance can cause errors. In view 
of this, on-line feedback control for the tracking problem, whereby the control decision is 
based on measurements of the current output error, appears more attractive. The technical 
literature on this subject is not very abundant. For tracking control, the Liapunov stability 
concept appears quite useful. Wang 8 used it to design a guidance law for a spacecraft docking 
with another spacecraft. The two docking objects are assumed to be three-dimensional rigid 
bodies and to have their own control system on board. Another assumption used in Ref. 8 
is that the motion of the target decays to zero with time. Recently, Novakovic® presented 
a technique using Liapunov-like methodology for robot tracking control problem. In this 
paper, the algorithm presented in Ref. 9 is adopted and modified for the tracking control of 
flexible space robots. 

In the case of flexible space structures, maneuvering motions excite vibration of the 
flexible members. There are two major control schemes for flexible manipulators. The first 
is based on linearized models derived from the nonlinear equations of motion of the flexible 
manipulator on the assumption that maneuver motions are much larger than elastic motions. 
Such a perturbation approach was developed by Meirovitch and Quinn 10,11 and applied by 
Meirovitch and Kwak 12,13 to the maneuvering and control of articulated flexible spacecraft 
and by Modi and Chang 14 and Meirovitch and Lim 15 to the maneuvering and control of 
flexible robots. The second is the adaptive control, 18 which does not need dynamical models. 
Instead, an auto-regressive-moving average (ARMA) model of system identification is used. 
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A control law for flexible manipulators based on the Liapunov method was proposed 
by Bang and Junkins. 17 It represents proportional and derivative control and includes a 
boundary force as a feedback force. This control scheme is valid only for problems in which 
the system approaches an equilibrium point in the state space. 

References 15 and 18 are concerned with flexible space robots of the type considered 
here, but the mission is more modest in scope. Indeed, in Ref. 15 the mission is to place a 
payload in a certain position and orientation in space and in Ref. 18 the objective is to dock 
with a target whose motion is known a priori. 

In this paper, a control scheme permitting a flexible space robot to track and dock 
with a moving target whose motion is not known a priori is presented. For the robot 
maneuver, on-line feedback tracking control is carried out by means of an algorithm based 
on Liapunov-like methodology and using on-line measurements of the target motion. For 
the vibration suppression, linear quadratic regulator (LQR) control in conjunction with 
disturbance compensation is carried out by means of sensor/actuator pairs dispersed along 
the flexible arms. A modified discrete-time control scheme is developed, and problems related 
to the digital implementation of the control algorithms are discussed. The control algorithms 
are demonstrated by means of a numerical example. 

2. EQUATIONS OF MOTION 

The flexible space robot and the coordinate systems are shown in Fig. 2. Body 0 repre- 
sents the robot base, assumed to be rigid. Bodies 1 and 2 are the robot manipulator arms 
attached in series to Body 0 and they are flexible. Body 3 is the end-effector/payload, also 
assumed to rigid. For planar motion, the robot base is capable of two translations, xo and 
yo, and one rotation, 0 q -, the two flexible arms are capable of the rotations 6 \ and &2 and the 
elastic vibrations ui and U2 and the end-effector is capable of the rotation #3. Referring to 
Fig. 2, the displacement vector Uo and velocity vector Vq for a typical point in Body 0 are 
as follows: 


Uo = R + Cq'Ro 

(3a) 

Vo = R + Cj'&oRo 

(36) 
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Similarly, for Body 1 


(4a) 

(46) 


Ui = R + C'o'Lo + C'i' (ri + ui) 

Vl = R + C'J’wqLo + Clu i (ri + ui) 4- Cf ui 


for Body 2 

U 2 = R + CffLo + Cf (Li + uu) + Cl (r 2 + u 2 ) 

V 2 = R + Cj’woLo + Cl 'uu (Li + u !2 ) + Cf uu 
+ Cfw 2 (r 2 4- u 2 ) + Cl u 2 

and for Body 3 

U 3 = R + clu + cl (Li + ui 2 ) 4- cl (L 2 + U 23 ) 4- Cl r 3 
V j = R+ Cj'woLo + Clu > i (Li + u 32 ) 4- Cfu 12 
+ Cfwa (Lj + u 23 ) + Cl u 23 4- Clu e r 3 


(5a) 

(56) 


(ба) 

( бб ) 


where 


fcos^ sin 0,1 i = 0, 1,2,3 

— sin 0* cos0iJ 

are matrices of direction cosines, 

r . ;l 

{ = 0, 1,2,3 


= 


0 -0i 

9i 0 


are 


skew symmetric angular velocity matrices, 

R=[x 0 yo] T , r 1 = [x 1 Of, r 2 = [x 2 of 


( 7 ) 

( 8 ) 
( 9 ) 


are position vectors and 

Ul = [0 Uif, u 2 = [0 u 2 f 
are elastic displacement vectors. Moreover, 

Ui 2 — Ui| tl = Li) u 23 U 2 |i 2 _Lj 


( 10 ) 


(ID 


5 



The elastic displacements are discretized as follows: 


Ui (x», t ) — (xj) ^ (t) , i — 1, 2 


(12) 


where $ j (x) are vectors of quasi-comparison functions 19 and f, (t) are vectors of general- 
ized displacements. Regarding the robot arms as beams in bending, the quasi- comparison 
functions c an be chosen as linear combination of the admissible functions 


, i AfcX A^x 

<pk = cosh — cos 


— <7* ^sinh — sin > A: = 1,2,... 


(13) 


L L 

which represent the eigenfunctions of a clamped-free beam for k odd and clamped-clamped 
beam for k even, where A* and are nondimensional parameters. 

Using Eqs. (3)-(13), the kinetic energy of the system can be written as 

T = / p.VfV.dD, = k r Mq (14) 

,=o ^ i= 0 •'Body » l 

'T 

where q = R r 6 q #i #2 #3 is the configuration vector and M is the mass matrix 

with entries given in Appendix A. 

The potential energy for the system is due entirely to the elasticity of the robot arms 
and can be written in the form 

2 


v = '£ = 2 V*q 


where 


K = block-diag [o K\ K 2 ] 


(15) 


(16) 


in which 

Ki = j* £/,$'' ($") T dxi, i = 1,2 (17) 

are the stiffness matrices for Bodies i, in which Eli denotes bending stiffnesses. Note that the 
gravitational potential is ignored here on the assumption that it represents a second-order 
effect. 

The control forces acting on the robot system include the horizontal and vertical thrusts 
F x and F y acting at the base center, the external torque Mq acting on the base, the internal 
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joint torques M\ y M 2 and M3 acting at the joints and the distributed internal moments T\ 
and r 2 generated by piezoelectric actuators on links 1 and 2. We define the control force 

ji 

vector as F = [i* 1 * F y Mq Mi M2 M3 rj r\ j . Then, the virtual work of the system 
can be written in the form 


8W =F x 8x 0 + FySyo + Mq89q + M\ ( 89\ — 89q) 

+ M 2 {892 ~ 89 \ - (Lx) 8Zx) + M3 ( 89 3 - 89 2 - ( L 2 ) 6£ 2 ) 

+ E *U*f (*«) *«, + E (*») = Q T *1 (18) 

»=1 1=1 

where Q is a generalized force vector defined as 


Q = GT 


(19) 


The entries of the matrix G are given in Appendix A. 

Lagrange’s equations for the system can be expressed in the symbolic vector form 

d_ (dT\ _ dT dV_ _ 

dt\dq) dq dq 


(20) 


Inserting Eqs. (15), (16) and (19) into Eq. (20), we obtain the system equations in the matrix 
form 

M(q)q + C(q y q)q + Kq= Q (21) 

The entries of the matrix C are also given in Appendix A. 

Equation (21) represents the equation governing the motion of the flexible space robot. It 
is used for computer simulation of the dynamical system. For the purpose of control design, 
Eq. (21) is conveniently separated into two sets of equations, rigid-body motion equations 
and elastic vibration equations. To this end, we write q = [q^ qjj and Q = [Q^ , 

rp rp 'T-'] X 

where q r = [zo yo #0 9\ 9 2 ^3] is a rigid-body displacement vector, q e = £2 j is 

an elastic displacement vector and Q r and Q e are corresponding generalized force vectors. 
Then Eq. (21) can be written in the partitioned matrix form 


M TT M re 
Mj, M, 


tt 


qr 
J Lq e J 


+ 


c rr c re 

G t r C et J 


qr 

Lq e 


+ 


0 0 

0 K 


h:]- 


Qr 

Q«J 


( 22 ) 
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After some algebraic manipulations, and ignoring higher-order terms in the elastic displace- 
ments, Eq. (22) can be separated into 

Mr far )qr + C T far , q r ) q r + d e fa, q, q) = Q r (23) 

and 

Me (qr) qe + C e fan qr) <3e + Ke far far far) qe + d r far far far) = Qe (24) 

where M r is the rigid-body part of the mass matrix M TT and C r is the rigid-body part of 
Crr- Moreover, M e = M ee , C t = C ee , K e consists of the stiffness matrix K and the part 
due to elasticity in the matrices M re and C Tt and d e and d r are disturbance vectors. The 
entries of the various matrices are given in Appendix B. The term d e in Eq. (23) is a linear 
combination of q e , q e and q e . It can be regarded as a disturbance due to the flexibility of 
the robot arms. The term d r in Eq. (24) is a function of q r , q r and q r . It can be regarded 
as a disturbance due to the rigid-body maneuvering of the robot. Equations (23) and (24) 
are coupled. The coupling between rigid-body motions and flexible vibration is provided 
in Eq. (24) by the persistent disturbance d r from the rigid-body motion, which causes the 
elastic motion q e , q e and q e . In turn, the elastic motion disturbs the rigid body motion 
through d e in Eq. (23). Equation (23) is used for the design of the maneuver control for 
tracking a moving target and Eq. (24) is used for design of control for vibration suppression. 

3. TRACKING CONTROL ALGORITHM USING LIAPUNOV-LIKE 
METHODOLOGY 

In this section, the general idea of Liapunov-like methodology for tracking control devel- 
oped for rigid robots 9 is introduced. 

The dynamical equation of a rigid robot i9 given by 

M (q) q + C (q, q) q = Q (25) 

and the kinematic relation between the robot configuration vector q and robot output vector 
y e is given by 

y.-f(q) (26) 
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so that 


y e = J (q) q 


(27) 


and 

y e = 7(q)q + j(q,q)q (28) 

where J (q) = [df/dq] is the Jacobian matrix. 

Because tracking is carried out by the end-effector, the tracking problem consists of 
driving the error e = yt — y e and its time derivative e to zero. To this end, a Liapunov 
function is defined by 

V = ^z r z (29a) 

where 

z = (e + /?e) (296) 

in which /? is a positive scalar. If the control is designed in such a way that 

V = — crV, (7 = In ^ /i 4 (30a, 6) 

where e is an arbitrarily small positive scalar and Vq is the initial value of V, it is guaranteed 
that the function V remains in the e-neighborhood of zero for t > t a , no matter how the 
target motion changes. This ensures that the error e and its derivative e are also very close 
to zero. 

We consider the nonlinear control law 


Q = M (q) u r + C (q, q) q 


(31) 


where u r is chosen in the form 


u r = w 


hi + hj 

z T Jw 


in which w is an arbitrarily chosen vector and 


(32) 


hi = z T (y t - jq + /3e) , h 2 = 0.5<rz r z = crV (33a, 6) 

It ran be shown that the control algorithm described above yields the desired result, i.e., 
Eqs. (30a, b). 
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The control algorithm possesses the following advantages: 

1) The control decision is made using on-line information of the current robot state (q, q) 
and target state (e,e and y<). The feedback control can automatically counteract ad- 
verse disturbances in space and achieve the final docking in an accurate and smooth 
way. 

2) The on-line calculation is relatively simple, as it involves neither inverse kinematics nor 
matrix inversions. 

3) Stability is always guaranteed by Liapunov stability theorem, as can be seen from 
Eqs. (30), no matter how the target motion changes. 

However, after applying the above algorithm directly to our space robot system and sim- 
ulating the system in both continuous time and discrete time, the results from discrete-time 
system exhibited some undesirable phenomenon, although the performance of the continuous 
system was good. As shown in Fig. 3, in which the solid line denotes continuous- time results 
and the dashed line denotes discrete-time results, the control force in discrete time exhibits 
periods of oscillatory behavior. Further numerical simulations show that the magnitude of 
the control force during chattering is bounded, although very large, and its mean value is 
close to the results of the corresponding continuous time system. Moreover, the occurence 
of the oscillating period is random, and the length of the oscillating periods and the length 
of the “good performance” periods are both unpredictable. This phenomenon is similar to 
the so-called “bursting”, which appears frequently in discrete-time adaptive systems and has 
been reported for almost a decade. 20 It is important to keep the control force from bursting. 
Otherwise the possibility exists that the control cannot be realized. To this end, a modified 
version of the above algorithm is presented, which also takes into account the flexibility of 
the robot arms. 

4. MODIFIED TRACKING CONTROL ALGORITHM FOR FLEXIBLE 
SPACE ROBOTS 

To apply Liapunov-like methodology to flexible space robots, we first extend the kine- 
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matical relation given by Eq. (26) to flexible space robots as follows: 

= 10-^0 sin 0 O + L\ cos 0\ + L 2 cos 0 2 + I3 cos £3 - u 12 sin 0\ - u 2 3 sin 0 2 
= tin 4- T.n mm 0n 4- T,t sin Qt 4- f/o sin 9i 4- L » sin 0% + ui<» cos 0t + Uia cos 0<> 


X e = Xo 

y e = yo + Lo cos 0o + L\ sin 0\ + L 2 sin 0 2 + L 3 sin #3 + ui 2 cos 0\ + u 2 3 cos 0 2 
0e = 0 3 

For kinematical analysis, we define q = 

rn — 

[“12 “23 
has the form 


( 34 ) 


q^ ql 


T* * - 

1 ] . The Jacobian matrix J, obtained 

r 

j = [Jr Ju\ 


Hr Huj » where qr was defined earlier and q tt = 
.ined by differentiating Eq. (34) with respect to q, 


(35) 


where 
J r = 


1 0 — Z-ocos^o — L\ sin 9\ — ui 2 cos 0\ — jL 2 sin 0 2 — u 2 3 cos 0 2 — Z- 3 sin 03 

0 1 — Losing Li cos^i — ui 2 sin^i L 2 cos 0 2 — u 2 3 sin 0% Z 3 CO 803 

- n 0 1 


0 0 


0 


Ju = 


Hence, we can write the relations 


— sin 9\ — sin 02 

cos cos 02 

0 0 

y e = Jk 

y e = + 


The dynamical equation for the rigid-body motion of the space 
We first define a nonlinear control law for Q r as follows: 


1 , 

(зба) 

(збб) 

(37) 

(38) 

robot is given by Eq. (23). 


Q r = M r (q r ) U r + C T (q r , qr) qr 


(39) 


Substituting Eq. (39) into Eq. (23), we obtain 

q r = u r - M- 1 d e (40) 

To prevent the bursting phenomenon, we propose a decoupled Liapunov function defined 

by 

Vi = Zi = €.{ + /?ej, i = 1,2,3 (41a, 6 ) 
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Taking the derivative of Eq. (41a) and using Eqs. (37), (38) and (40), we obtain 

V% = Z{hi — z{ ([</rU r ]j ~ *d e J^ > * = 1,2,3 (42) 

where [ ]j denotes the i — th element of a vector and hi are the components of the vector 

h = y* — jq + /?e — J tt q tt (43) 

Because M r is a positive definite matrix, M~ l is bounded, and we note that J r is also 
bounded. Moreover, from Eq. (B.3) in Appendix B, we see that d e is a linear combination 
of q e , q e and q e . We then assume that d e is bounded in accordance with our ultimate goal 
of vibration suppression. Hence, we can assume that the term J r Af~ 1 d e j . is bounded and 
satisfies the relation 

. = 1,2,3 (44) 

Prom Eq. (44), we have 

zi [j r M r _1 d e ] . < t = 1,2,3 (45) 

If we can determine a vector u r that satisfies the following conditions: 

Zi [J r u] r = zihi + + |z;|6i, t = 1,2,3 (46) 

then 

Vi=^atiz? +\j r M- 1 d e .-\zi\6i < -^otiZi = -aM, . = 1,2,3 (47) 

According to the Liapunov stability theorem, Eq. (46) is the sufficient condition for our 
tracking problem. We further simplify Eq. (46) by assuming Zi ^ 0, so that 

[J r u r ], = h t + ^atiZi + sgn(z,)^, i = 1,2,3 (48) 

or 

[J r u r ]j = 3 t , . = 1,2,3 (49) 

with 

Si = yu ~ [-M] . + foi ~ [«/«q«L + + sgn (- 2 .) (50) 
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Equation (49) can be expressed in the matrix form 

J r u r = s (51) 

where s = [si sj 33]^ and J T is a 3 x 6 matrix. The solution of Eq. (51) does not yield 
a unique u r . This agrees with Eq. (32) in the original control scheme in which w is an 
arbitrarily chosen vector. Here we can simply prescribe the redundant degrees of freedom 
and then solve Eq. (51) accordingly. 

As a simple example, we constrain three components of u r by taking 

^r3 = tir4 = ^r5 — 0 (52) 

for the entire tracking period and use Eqs. (51) to solve for the other three components of 
u r on-line, with the result 

U r i = + L3Sin#3U r 6 

Ur2 = 32 - L 3 COS 0 3 Ur6 (53) 

Ur6 = -S3 

The above algorithm for u r , together with Eq. (39), represents the maneuver control for 
a flexible space robot tracking a moving target whose motion is not known a priori. The 
control algorithm requires that the following conditions be satisfied: 

1) The output error vector e and its time derivative e can be measured on-line. 

2) The target output acceleration y* can be measured or estimated on-line. 

3) The robot rigid-body displacement vector q r and its time derivative q r can be measured 
on-line. 

4) The elastic tip displacement vector q u and its time derivatives q tt and q u can be measured 
on-line. 

5) The elastic vibration of the robot axms should be controlled so that a reasonable value 
for the upper bound can be set. 

In addition to the advantages of the original algorithm mentioned in Sec. 3, the modified 
control algorithm presented here provides two extensions from the original one. 9 The first ex- 
tension is that the flexible effect of the robot arms is incorporated into the control algorithm. 
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It is reflected in the kinematic relations expressed by Eqs. (34) and in the term sgn (z») Si in 
Eq. (50) which is associated with the vibration disturbance vector d e in Eq. (23). The second 
extension consists of the use of decoupled Liapunov functions, Eqs. (41), to eliminate the 
bursting phenomenon (Sec. 3) when the control algorithm is implemented in discrete-time. 

5. VIBRATION CONTROL 

Because of coupling between the rigid-body motions and the elastic vibration, the per- 
formance of the tracking control depends on how well the vibration suppression is carried 
out. Without vibration control, the tracking cannot be truly realized for a flexible space 
robot. Our objective is to drive the elastic motion state q e , q e close to zero during the 
tracking and docking operation. We recall that the motion of the elastic vibration of the 
space robot is described by Eq. (24), which represents a linear time- varying system with a 
persistent disturbance term d r due to the rigid-body motions. 

We propose to control the vibration in discrete time. To this end, we separate the 


generalized control force Q e into 

Q e (k) = Q er (k) + Q„(k) (54) 

The discrete-time control algorithm for disturbance compensation is expressed by 

Qer (k) = d r (q r ( k ) , q r (fc) ,q r (fc)) (55) 

If the disturbance is cancelled perfectly, Eq. (24) becomes 

M e (q r ) q e + C e {q r ,q r ) q e + K e (q r ,q r ,q r ) q e = Q ee (56) 

[ rp JH T' 

q e (k) q e (k) be the state vector and u(fc) = Q ee (k) the control 
vector, the discrete-time state space counterpart of Eq. (56) can be written as 

x (k + 1) = A (k) x (k) + B (k) u (k) (57) 

where the coefficient matrices are given by 

A ( k ) = e A ( kT \ B (Jfe) = (e A{k7 "> - /) A T (kT) B (kT) (58a, b) 
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in which 


-M7 l K e -M~ l C e . 


The performance index for the discrete-time LQR is given by 21 


J = \ £ [x T (*) Q (k) + u r (fc) Ru (fc)] 


yielding the control law 


u(Jfc) = - (R + B(k)K(k)B(kj) 1 B T (k)k(k)A(k)x(Jc) 


(59a, 6) 


(60) 


(61) 


where K (k) satisfies the discrete-time algebraic Riccati equation 

K (k) = A t (k) [k (k) - K (k) B (k) (r + B t (k) K (t) B (*)) B T (k) K (k) ] A (k) + Q 

(62) 

Direct application of the discrete- time control algorithm described by Eqs. (55) and (61) to 
our problem causes severe instability. The reason is that the discrete-time control force Q er 
in Eq. (55) is not able to cancel the continuous disturbance term d r in Eq. (24) perfectly. 
Hence, the LQR control design based on Eq. (56), in which the disturbance is absent, is no 
longer appropriate. The error accumulates with time and it finally results in instability. To 
resolve this problem, a modified vibration control algorithm is proposed in the next section. 

0. MODIFIED DISCRETE- TIME VIBRATION CONTROL ALGORITHM 

An examination of the disturbance term d r in Eq. (B.14) of Appendix B, i.e., an exam- 
ination of 

d r = q r + C er q r (63) 

reveals that q,. in the first term is the major cause of the system instability. Usually q r ( Jc ) 
is not available and q r (& — 1) is used as an estimate of q r (fc). Stable performance of the 
system can be achieved only if q r (&) can be measured or estimated perfectly. Even a very 
small error in q r appearing in Eq. (63) can result in failure of the LQR design. To avoid 
use of q r in Eq. (63), we replace q r by u r , so that the disturbance compensation scheme 
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becomes 


Qer (k) — d r (q r (fc) , C[r (&) , U r (fc)) 

(q, (ft)) u r ( k ) + C„ (q r (k) , q. (*)) q, (k) (64) 

where u r (k) is calculated by the tracking control algorithm given by Eq. (51). We then 
substitute Eqs. (63), (64) and (40) into Eq. (24) and obtain the system equation as follows: 

Me(qr)q e + C' e (q r ,q r )q e + (q r , q r , q r ) q e - M^M r _1 d e = Q ee (65) 

As shown in Appendix B, d e can be expressed as 

d e = Mrt q e + C re q e + (Klf + Kh) q e (66) 

where K e M and Kq are given by Eqs. (B.6) and (B.8), respectively. Substituting Eq. (66) 
into Eq. (65), we obtain the modified linear time- varying system 

M* (q r ) q e + C* (q r , q r ) q e + K* (q r ,q r ,qr) q e = Q« (67) 

where, comparing Eqs. (56) and (67), we observe that matrices M*, Cl and A - * represent 
modified coefficient matrices given by 

Ml = M e - Mj e M~ l M re (68a) 

Cl = C 6 - Mj e M~ l C rt (686) 

Kl = K e - ( K e M + K e c ) (68c) 


Based on Eqs. (67) and (68), we can follow the same procedure as in Sec. 5 and obtain the 
control law for Q M . The simulation results using the modified control scheme showed stable 
performance. F\irther numerical simulations showed that even in the case of a system with 
only the mass matrix M e modified, i.e., a system described by 

Ml (q r )q« + C e (q r , q r ) q« + K e (q*,q r ,q,) q« = Q« (69) 

the LQR control law is still able to produce good system performance. This is because the 
first term on the right side of Eq. (66) is dominant, so that using C e and K e instead of Cl 
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and K*> respectively, is equivalent to dropping the second and third terms in Eq. (66), which 
does not affect the system performance very much. 

7. NUMERICAL EXAMPLE 

We assume that the parameters for the flexible space robot shown in Fig. 1 have the 
values 

mo = 40.0 kg, mi = m2 = 10.0 kg, m3 = 2.0 kg 
Lq = 2.5 m, L\ = L2 = 10.0 m, L3 = 2.0 m 

(70) 

S t = Sy = 0, I s = 83.333 kg m 2 , /„ = 333.333 kg m 2 
Eli = EI 2 = 10 4 kg m 2 


The target motion is not known a priori and must be measured on-line. However, for 
simulation purposes, we choose an example target trajectory as follows: 

x t (t) =10.0 sin 

y t (t) =10.0 -I- 10.0 sin , t 6 [0,5.0 s] (71) 


The initial conditions of the space robot are given by: 

q f (0) = [0-0 - 15.0 0.0 0.5 tt 0.4775a- 0.25jr] r , q r (0) = 0 
q e (0) = [0.01 ... O.Olf, q e (0) = O 
The parameters of the control synthesis design are 


(72) 


0 = 20.0, £=10" 3 , i, = 2.5 j, 0i=2O, i = 1,2,3 


(73) 


We designate the three redundant degrees of freedom in u r as u f 3, u r 4 and u r s. They 
are defined for two different cases as follows: 

Case 1: 

U r 3 = ^r4 = ^rS ~ 0, 
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0 < t < 5 s 


(74) 



Case 2: 


0, 

t < 0 

4A0 O /t}, 

0 < t < tf/2 

-4A0oA/> 

tfl 2 <t<tf 

0, 

t > t f 

0, 

t < 0 

4A exit), 

0 < t < tf/2 

-4A0i /*}, 

tf/2,t<t f 

0, 

t > tf 

0, 

t < 0 

4A0 2 

0 < t < tf/2 

— 4A e 2 /t), 

tf/2 <t<t f 

0, 

t > t f 



(75a) 


(756) 


(75c) 


where tf = 4.0 s, AOq = | rad, A6\ = j rad, and A6 2 = — f rad. 

For a rigid space robot, Eqs. (74) and (75) represent constraints on the acceleration of 
the robot configuration. In Case 1, the mission amounts to keeping the base attitude 


and the two joint angles 6\ and 9 2 constant while tracking a moving target. In Case 2, 


the mission implies bang- bang maneuvers involving a base attitude change of A0 q and arms 


angle changes of A#i and A9 2 while tracking a moving target. 


The constraints cannot be realized perfectly for a flexible space robot due to disturbance- 
causing vibration. However, the performance can be improved by vibration control. Because 
the major objective here is to track the moving target, we use the constraint equations, 


Eqs. (74) and (75), to eliminate the robot redundancy. 

For vibration control, the LQR design parameters are chosen as 


R = diag [/ n xn /nxn 


(76) 


Q = diag [2.0 x 10 4 / nX « 10 4 / nxn 2.0 x 10 4 / nxn 10 4 / n 


xn 


The elastic displacement for each of the two arms was modeled by means of five quasi- 
comparison functions. 19 

The system performance under the tracking and docking maneuver is simulated over 5 s. 
To this end, the tracking control algorithm presented in Sec. 4 and the vibration control 
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algorithm presented in Sec. 6 axe used. The simulation is performed in discrete-time with a 
sampling period T = 0.001 s. 

Figures 4a and 4b show time-lapse pictures of the robot configuration for Cases 1 and 2, 
respectively. For Case 2, time histories of the tracking error e and its time derivative e axe 
shown in Figs. 5a-5c, time histories of the tip elastic displacements of the two flexible links 
are shown in Fig. 6 and time histories of the control forces and torques for the rigid-body 
maneuver axe displayed in Figs. 7a-7c. Time histories of the control torques acting on the 
flexible bodies for disturbance rejection and LQR control are shown in Fig. 8 and Fig. 9, 
respectively. The results axe very satisfactory, with control achieved in less than one second. 

8. SUMMARY AND CONCLUSIONS 

This paper is concerned with the control of a flexible space robot executing a docking 
maneuver with a target whose motion is not known a priori. The control is based on on- 
line measurements of the target motion. The dynamical equations of the space robot axe 
first derived by means of Lagrange’s equations and then sepaxated into two coupled sets 
of equations suitable for rigid-body maneuvers and vibration suppression. Controls for the 
rigid-body maneuver and vibration suppression are developed and implemented in discrete 
time. Problems arising from digital implementation of the control algorithms axe discussed. 
Then, modifications of the control algorithms so as to prevent the problems axe made. 

The control scheme presented can be applied to two-dimensional, as well as three- 
dimensional problems. Furthermore, it has the flexibility of solving different problems by 
defining appropriate output vectors other than the end-effector output vector. For example, 
if the mission involves tracking and docking with an orbiting target while its base attitude is 

rp 

to be kept constant, we can define the output vector as y e = [ x t y e 9 t 0o] and the target 
output vector as y ( = [x t y t 6 t 0] T , and then the proposed tracking control algorithm can 
be used to drive the error vector e = y t — y e and its time derivative e to zero. 

A numerical example is used to demonstrate the control scheme. The simulation results 
have shown very good system performance in both the tracking maneuver and the vibration 
suppression. 
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APPENDIX A - Matrices in the Equations of Motion 

The mass matrix M appearing in Eq. (14), as well as in Eq. (21), is defined as 
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in which 


(A.l) 


(A.2) 


oi = -Stisi - *S*iCi, a 2 = -S t2 S2 - $f 2 $ 2 c 2 

03 = ‘S'tici - ^n£i s i> a 4 = S t2 c 2 - $a£ 2 s 2 

as = + *S*i£ocio, = St 2 Los 2 o + $^£2^0020 

a 7 = St 2 L\c 2 \ + S t 2Qi 2 Zi32i - $at 2 Lis 2 i + ^S^ 2 $ ? 2 ^i c 2 i (A.3) 

rry m 

as = S3L1C31 + 53$12^1 5 31, a 9 = ^3^2^32 + ^3^23^2^32 

bi = $t2$lUl 5 21> b 2 = -*i2#S$2 a 21 

— m • rjr\ 

h = hi + £1 m77^i, h = hi + i 2 m 88$2 
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and 

mn = — $tisi, m27 = *tici, ITI 37 = $ti£osio 

ITU7 = 4^1 4- (m2 + m3) Z/1$12, ITI57 = St 2 $UC 21 , m 67 = «S3$12C 3 1 

mi 8 = -$<232, IT128 = $<2C2, HI 3 8 = $12C21, ni67 = 5 , 3$12C31 

(A.4) 

rri 48 = $<2^1021, ni58 = $2 + m3.Z/2$23, ni68 = ■S , 3$23C32 
m77 = Ai + (m2 + m3) $12$12> m 78 = $12$«C21 

'P 

mg 8 = A2 + m3$23$23 

and we note that Sj = sin 0<, c\ = cos 6 { , sjy = sin ( 9 { — 0y) and c^y = cos (0< — 6 j). Moreover, 
we have used the following definitions: 


mj = mo + mi + m2 + m3 

Stz = So* sin 0o + So y cos 0o + (mi + m2 + m3) Lq cos 0 o 
Sty = — •S'oz COS 00 + Soy sin 0 o + (mi + m2 + m3) Lq sin 0 o 

Sti = S\ + (m 2 + m3)Li, S t 2 = S 2 + mzLi 
ho — hz + loy + (mi 4- m2 + m3) Lq 

0 o 

hi = h + ( m 2 + m3) L\ y It2 = / 2 + m3L 2 

$tl = $1 + ( m 2 + m 3 ) $12, $t2 = $2 + m 3 $23 

in which 


m = / B0 dy i PiiDi i : 

= 0,1, 2,3 

Si = /Body , PiXidDi ' 

U = Liy^ iD " ‘ 

$ 0x Jb ody 0 

PqxcIDq , 
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'“‘“•feodyo'’ 01 ’^ 0 ' 
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./Body • ’ 
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$23 = $2 (X 2 ) | X2=ij 
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The matrix G in Eq. (19) is defined as 

'1 0 0 0 0 0 0 T 0 T ' 

0 1 0 0 0 0 O r O r 

0 0 1 -1 0 0 O r O r 

^ 0 0 0 1 -1 0 o r 0 T 

G ~ o 0 0 0 1 -1 o T o T 

0 0 0 0 0 1 o r o r 

0 0 0 0 0 G\ 0 

.0 0 0 0 0 -$' 2 (L 2 ) 0 g 2 . 

where primes denote spatial derivatives and 

[»f(x a ) ... •?(*«»)] * = 1.2 


(A.7) 


(A.8) 


in which m is the number of actuators on each link. Here m is equal to the number of modes 
and Gi are square matrices. 

The coefficient matrix C in Eq. (21) is defined as 

'0 0 C13 C14 Cis Ci6 C17 Cia 

0 0 C23 C24 C25 C 2 6 C 2 7 C*28 

0 0 0 C34 C35 C36 C37 C38 

p _ 0 0 C43 0 C45 C46 C47 C48 f A.9) 

0 0 C53 C54 0 C56 C57 C58 

0 0 C63 C64 C*65 0 C67 ^68 

0 0 C73 C74 C75 C76 0 C78 

.0 0 C83 C84 Css Cb6 C *87 0 . 

where 


Gn = StyOo, Cu = (-■S'tici + b\ } C 15 = (-5 t 2C2 + h 

C16 = — 830363, = — 2 $jjCi 6 i i, Cj 8 = — 2 $^C 2^2 

C 2 3 = -Stzbof Cu = (—■S'tl^l - $tl£l c l) ^1 5 ^25 = {—St 2^2 - ^a^2 c 2) &2 

C 2 6 = —S 3 S 3 O 3 , C 2 7 = — 2$4jSi0i, C28 = — 2$^32^2 

C34 = (-S'tiXocio - * 5 ei^«io) bi,C 3 s = {St 2 ^oc 2 o - 02 

C36 — S3LQC3063, C37 = 2 $^Locio^i> C38 = 2 . 0 $r 2 L 0 C 20^2 

C43 = (-SnLocio + $ti£i-kosio) 0o 

C 45 = {—St 2 L\S2\ — 2 L\C21 + 5t2^!2^1 c 21 — ^S^2^12^1' s 2l) 02 
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(746 = (-SiLiSil + S , 3 $i 2 ^ 1 c 3 l) ^3 
C 47 = 2tf (Ai + (m 2 + m 3 ) *i2*f 2 ) *1 
C48 = 2 (“Ll321^S + *12£l c 21*&) ^2 

C53 = (— St2L<)C2Q + ^^2^0-S2o) #0 

Cii = (5(2^1321 + ^«^2^1 c 21 - St2»fa*l<*l + #1 

C56 = (-S3L2S32 + *S'3^M$2 c 32) # 3 . C57 = 2 ( 5 t 2321^!2 + ^tt&^l^ia) ^1 

C58 = 2^ (^2 + m 3 $23$r 3 ) ^1> (^63 = -S3L0C30O0 

(^64 = (S 3 L 1 S 31 ~ ■S , 3^l2$l c 3l) #1. Ces = [S 3 L 2332 - 5 3 $23^2 c 32) #2 

C67 = 25' 3 33i$23^2) ^68 = 253332$ 23^2 > C73 = — ^tl-^0 c 10^0 

C74 = - (Ai + (m 2 + m3) #„•£) > C75 = (-■S'«321^12 _ $t2$2 c 21$12 
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Cs4 = (la3 2 l$t2 ~ $l2£l c 21^t2) #1> ^85 = ~ (A2 + m3$23$23) ^2^2 

Cs6 = —S 3332 & 2303 , Cfo = 2$t2$12’ s 2101 


APPENDIX B - Matrices in the Partitioned Equations of Motion 


The mass matrix M r and the coefficient matrix C r in Eq. (23) are defined as 
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The disturbance vector d e in Eq. (23) is defined as 
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and the coefficient matrix K t is defined as 


K e = K + Km + K c 
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Fig. 3 - The Bursting Phenomenon 







Fig. 4a - Time-Lapse Picture of the Robot Configuration - Case 1 



Fig. 4b - Time-Lapse Picture of the Robot Configuration - Case 2 
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Fig. 5b - Time History of the y-Component of the Tracking Error and Tracking Error 
Rate - Case 2 
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Fig. 7b - Time History of the Control Torques for the Rigid-Body Rotation of the Robot 
Base and Body 1 





Fig. 7c - Time History of the Control Torques for the Rigid-Body Rotation of Bodies 2 
and 3 
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Fig. 8 - Time History of the Control Torques Acting on the Flexible Bodies for Dis- 
turbance Rejection 
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Abstract 

A variety of engineering systems, such as automo- 
biles, aircraft, rotorcraft, robots, spacecraft, etc., can 
be modeled as flexible multibody systems. The individ- 
ual flexible bodies are in general characterized by dis- 
tributed parameters. In most earlier investigations they 
were approximated by some spatial discretization pro- 
cedure, such as the classical Rayleigh-Ritz method or 
the finite element method. This paper presents a math- 
ematical formulation for distributed-parameter multi- 
body systems consisting of a set of hybrid (ordinary 
and partial) differential equations of motion in terms 
of quasi- coordinates. Moreover, the equations for the 
elastic motions include rotatory inertia and shear de- 
formation effects. The hybrid set is cast in state form, 
thus making it suitable for control design. 

1. Introduction 

A problem of current interest is the dynamics and 
control of multibody systems. Indeed, a variety of en- 
gineering systems, such as automobiles, aircraft, ro- 
torcraft, robots, spacecraft, etc., cam be modeled as 
multibodies. In many engineering applications the bod- 
ies can be assumed to be rigid (Refs. 1-12). In many 
other applications, the flexibility effects have to be in- 
cluded (Refs. 13-24). For the most part, flexible bodies 
have distributed mass and stiffness properties, which 
is likely to cause difficulties in producing a solution. 
As a result, it is common practice to approximate dis- 
tributed systems by discrete ones through spatial dis- 
cretization, which cam be carried out by means of the 
classical Rayleigh-Ritz method or the finite element 
method (Ref. 25). The discretization process amounts 
to elimination of the spatial coordinates. The equations 
of motion for the discretized system are derived quite 
often by the standard Lagrangian approach. For more 
complex motions, an approach using quasi-coordinates 
seems to offer many advantages (Refs. 26-29). 

1 Supported by the AFOSR Research Grant F49620- 
89-00045 monitored by Spencer T. Wu and by the 
NASA Research Grant NAG-1-225 monitored by Ray- 
mond C. Montgomery. 

* University Distinguished Professor. Fellow AIAA. 
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Quite recently, there has been some interest in 
working with distributed models as much as possible, 
thus avoiding truncation problems arising from spa- 
tial discretization. Consistent with this, hybrid (ordi- 
nary and partial) differential equations of motion have 
been derived for flexible multibody systems in Refs. 30 
and 31, using the approach of Ref. 25. Hybrid equa- 
tions of motion in terms of quasi-coordinates have been 
derived for the first time in Ref. 26 for a spinning rigid 
body with flexible appendages and generalized later in 
Ref. 32 for a flexible body undergoing rigid-body and 
elastic motions. This paper extends the general theory 
developed in Ref. 32 to systems of flexible multibod- 
ies. In addition, the equations for the elastic motions 
include rotatory inertia and shear deformation effects. 

2. Kinematics 

We are concerned with structures consisting of a 
chain of articulated bodies t (i = 1,2 ,...,#), which 
implies that two adjacent bodies i — 1 and t are hinged 
at Oi (Fig. 1). To describe the motion of the system, 
it will prove convenient to conceive of a set of body 
axes ZiyiZi with the origin at 0* and attached to body 
i in undeformed state. The bodies are assumed to be 
slender, with axis z* coinciding with the long axis of 
the body. As the body deforms, z* remains tangent to 
the body at Oi. At the same time, we consider another 
set of body axes referred to as intermediate axes, 

with the origin at 0% and attached to body t — 1 so that 
z[ is along the long axis. We will also find it convenient 
to introduce an inertial frame of reference XY Z with 
the origin at O. 

We denote the position vector of point Oi relar 
tive to the origin O by Roi = [X oi Y ai Z ot ] T . Then, we 
denote the position of a typical point Vi in the unde- 
formed i body relative to O; by r< and the elastic dis- 
placement oiVi by u». Hence, the radius vector from 
O to Vi in displaced position is simply 

R< = C'Roi + ti + Uj, i = l,2,...,N (1) 

where C* is the matrix of direction cosines of axes 

with respect to axes ij/»_ i «»— 1 , and note that 
the vector is in terms of components along the body 
axes and the vectors R*, r, and u< are in 

terms of components along the body axes ZiS/iZ,. 
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We consider here bodies in the form of bars with 
the long axis x, parsing through 0{ and 0,+ i when 
the bara are undeformed.' We are concerned with bar® 
undergoing torsion about axil x* and bending about 
axes Vi and Z*, u well ai shearing distortion in the 
and Zi directions. Then, the vectors r, and u« can be 
written in the more explicit form 

r< = [*i 0 Of, Ui(*i,t) = [0 uy i(x it t) 

(2a,b) 

The radius vector R<* depends on the motion of the 
preceding t — 1 bodies in the chain. In particular, we 
can write the following recursive relation: 

R« — C^i R-o,»- 1 "hr$- i(A- i ) ~b i , t), 

* = 2,3, ...» W (3) 

where is the length of body i — 1. Note that 
R 0l = R 0 i(t) is simply the radius vector from O to 
the origin O i of the body axes of the first body in the 
chain. 

At this point, we propose to define the rotational 
motions. In the first place, it will prove convenient to 
introduce a set of body axes £*^6 attached to a typical 
beam cross section originally in the nominal position Zi 
and moving with the cross section as body i deforms. 
In this regard, note that 

coincide with x<y-z<. Then, denoting the angle of twist 
by and the bending rotation angles by and t 
we conclude that axes experience the angular 

displacement 

= [*/’«(*», t) ( 4 ) 

with respect to axes z^z,. On various occasions 
throughout this paper, we encounter skew symmetric 
matrices derived from vectors. As an example, if a 
typical vector r has components x, y and z, then the 
associated skew symmetric matrix has the form 


Then, the matrix of direction cosines of axes x^x* 
relative to axes is simply 

Q = CiEi-iiii-ut) (7) 

From kinematics, the velocity vector of the typical 
point Vi in displaced position in terms of the rotating 
body axes z^z^, has the expression 

V* =V* + + u i) + Vi 

=V« + (r i + Ui) T n ri + v i , i = 1 , 2 , ... t N ( 8 ) 

where is the velocity vector of the origin 0{ y 
is the angular velocity vector of axes ZiyiZi relative to 
axes XYZ and 

= Ut (*t,0 (9) 

is the elastic velocity vector relative to z^z,, all in 
terms of z^z^ components. We note that the velocity 
vector of point 0* can be written in the recursive form 

V* sQVi-iCA-i.t) 

-Cl {v 0li -i + [r«-i(4-i) + fii-i( 4 -i,t)] r Or,*- i 

+ v<_!( A-i.O). < = 2 , 3 ,. .., 1 V ( 10 ) 

Moreover, introducing the notation 

fi«t(x<,t) = i = >N (11) 

the angular velocity vector of the cross-section ad axes 
&>7iCi relative to the inertial apace is simply 

n< = n r< + n*(x<,t), < = 1 . 2 , . . . , n ( 12 ) 

Finally, letting u>i be the angular velocity vector of 
axes z^iZi relative to axes x^y-z-, in terms of r,y,Zi 
components, the angular velocity vector of is 

given by the recursive formula 


0 


r 


z 



-y z 0 


( 5 ) 


In view of this definition, the matrix of direction cosines 
of relative to z^z* can be shown to have the 

expression 

Ei(*i,t) = / - (6) 

in which I is the 3x3 identity matrix, and we note 
that Eq. (6) follows from the assumption that the 
components of *** ®mall. Next, we assume that 
axes z^iZi are obtained from axes x'y^x- through the 
rotations where j can take the values 1, or 1,2, 
or 1,2,3, depending on the nature of the hinge at Oi 
and denote by (7,(0,) the matrix of direction cosines 
of z^iZi relative to zjy^z-, where Si = [$u 9^ 9iz] T . 


f^ri — Cl — 1 (^t- 1 1 1) + Wi 

= 1 “b 1 (f*- 1 j 0] "b W,, 

1 = 2,3,. ...AT (13) 

where the second equality follows from Eq. (12). 

3. Standard Lagrange's Equations for Flexible 
Multibody Systems 

The motion of our multibody system is described 
in terms of rigid-body displacements of sets of body 
axes and elastic displacements relative to these body 
axes. As a result, the equations of motion are hybrid, 
in the sense that they consist of ordinary differential 
equations for the rigid-body displacements and partial 
differential equations for the elastic displacements. The 
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equations of motion can be derived by means of the 
extended Hamiltion’s principle (Ref. 33), which can be 
stated in the form 


L 


(SL + 6W)dt = 0, Sq = 0, <Su, = Sj { = 0, 

1 = 1,2,...,// &t t = ti,t 3 (14) 


where 


L = T-V 


(15) 


is the Lagrangian, in which T i s th e kinetic energy and 
V is the potential energy, and SW is the virtual work. 
Moreover, q is the rigid-body displacement vector, and 
u*, (* = 1,2,...,//) are the elastic displacement 

vectors introduced earlier. Hence, before we can derive 
equations of motio n, we must derive general expressions 
for T, V and SW. 

Taking the x»-axis to coincide with the centroidal 
axis of the undeformed beam, the kinetic energy can be 
shown to consist of two parts, one due to translations 
and one due to rotations (Ref. 25). Hence, using Eqs. 
(8) and (12), the kinetic energy can be expressed in the 
form 

* rli 

T = Y, / Tidxi 

where 


(16) 


t = l -( P iV7Vi + nf SdOi) = \\piV T oi Voi 

-T 

+ Jiflri + P % Ui + flri 

+ 2p<V£u* + 2fl^ t 5*Ui + J ci^ri 

+ + 2n£/ciV\] 

=\\fiiVlv oi + n£/ t< n r< 

— T 

+ Pilljui + Jd^i + 

+ 2p i VZii i + 2n? i (hiXi + jM} (17) 

is the kinetic energy density of member i, in which p { 
is the mass density and 


Jti — J\ Jc i 


(18) 


is the total moment of inertia density matrix, where 


Ji =Pi(n + + Ui) 


r..\T 


r «2i + «5» -*<%< 


-Pi 


v* 

— XjUyi 


®i + 


Xj Uji 
“ Uyitijj 

+ U £i 


and 


= diag(J, 


'v»»» 




(19a) 


(19b) 


in which JdMii and Jm%m% are cross-sectional mass 
moments of inertia densities, and note that, because the 
elastic deformations are relatively small, they are ap- 
proximately equal to and respectively. 

Moreover, Si is obtained from 

S« = pt(rj + Ui) — p\\%i Uy,' tij,] (20) 

which is recognized as the first moments of inertia 
density vector. 

Assuming that differential gravity effects are neg- 
ligibly small, the potential energy reduces to the strain 
energy. As indicated earlier, the elastic members un- 
dergo torsion about Xi and bending about y* and z», 
as well as shearing distortions in the y* and z* direc- 
tions. Referring to Fig. 2, we conclude that the rela- 
tions between the bending displacements Uyi and u Mii 
the bending angular displacements rj> y i and and the 
shearing distortion angles f} y i and are 

u y» = u j« = “Vv* “ Pyi (21a»b) 

where primes denote partial derivatives with respect to 
x*. From mechanics of materials, the relation between 
the twisting moment M ti and the twist angle rp ti is 
simply 

M gi = k zi GiI zi ^ mi (22) 

where kd is a factor depending on the shape of the cross 
section and Gild is the torsional rigidity, in which Gi is 
the shear modulus and Id is the polar area moment of 
inertia about axis x*. Moreover, the bending moments 
are related to the bending rotational displacements by 

M y i = Eilyrfyi, M Mi = EiI ti tV Mi (23a, b) 

in which Ei ia Young’s modulus and lyi and I,i are 
area moments of inertia about axes parallel to y* and 
Zi, respectively, and passing through the center of the 
cross-sectional area, and the shearing forces are related 
to the shearing distortion angles according to 

Qy% — kgiGiAiftgi , Q M i = —kyiGiAiflyi (24a, b) 

where kyi and k M i are factors depending on the shape 
of the cross sectional area, Gi is the shear modulus and 
Ai is the cross-sectional area. 

The strain energy can be expressed as 



where, using Eqs. (21)-(24), 

Vi =- (AfffiV’ci + S^yilfiyi + A + QyiPti "* QsiPyi) 

= \ [k'iGilMi) 7 + EJyWytf + EilMi) 2 
+ k yi G i A i (u l yi - if+k'iGiAiMi + t/. yi ) J ] (26) 


3 



i s the potential energy density for member i. 

Next, we wish to develop an expression for the 
virtual work due to nonconservative actuator forces and 
torques. Using the analogy with Eqs. (8) and (12), the 
virtual work can be written in the form 


N r & 

6W = +mffi©-)dxi 


»=i 


N 


+ '£m'J691 


i = 2 


=E{/V( 5R ^+ f i’ fi0 *‘+^) 

i N 

+uf(6&; i +6il> i )]dxi \ + J2 M '* Se ‘ 

* 1 = 2 

= f;[F;, r 5RV+M;f50V 

1=1 <■ 

r ti 1 JL 

+ 1 (f? 6\ii+m[ tfV’i) dxi +^M;, T (50* (27) 

J 0 1 = 2 


in which f* and mi are distributed actuator forces and 
torques acting over the domain i, are torque ac- 
tuators located at points and acting on both mem- 
bers s - 1 and i, for t = 2, 3, . . . , N , <5R- is the virtual 
displacement vector of point V{ y <50- is the virtual ro- 
tation vector of axes <50* is the virtual rotation 

vector of axes relative to axes r'y-z-, <5RV is 

the virtual displacement vector of point O t and <50^ 
is the virtual rotation vector of axes z»y*Zt relative to 
axes XYZ ) where all of these vectors are in terms of 
components along axes ZiyiZi, and asterisks indicate 
quasi-coordinates (Ref. 33) and associated forces and 
torques. Note that the term f J uj 8®* ri was omitted 
from £RJ on the basis that it is second-order in mag- 
nitude. Moreover, 


rli rU 

F' ri = / f idx it = / (fit + mi) dXi (28a, b) 

Jo Jo 

are, respectively, resultant forces and torques acting on 
member i. 


respectively, are also all independent. It will prove con- 
venient to introduce the rigid-body motion vector 

q(0 = [RfiW (0 (*)••• *£(‘)] T (29) 

so that we propose to derive a vector Lagrange ordi- 
nary differential equation for q(f ) and N pairs of vector 
Lagrange partial differential equations for u*(z*,t) and 
VnO^iO (* = 1,2,..., N ). To this end, we wish to ex- 
press the Lagrangian in general functional form, and we 
note that the Lagrangian contains not only q, u* and V>» 
but also time and spatial derivatives of these vectors. 
Moreover, we observe from Eqs. (3), (7), (10) and (13) 
that the Lagrangian contains terms involving u »(£,*)» 
Ui(4,<), V\(4,t) and V\(4,*). Such terms will con- 
tribute to the dynamic boundary conditions accompar 
nying the partial differential equations for u*(xi,i) and 
, t). In view of this, we express the Lagrangian in 
the general form 

L =£-(q,q ) u <> u-,u< l V»i.V><.V><.u*(A,0. 

Ui(4,*)>V\(4,t).V\(4. t )] ( 30 ) 


The extended Hamilton’s principle, Eq. (14), calls 
for the variation of the Lagrangian, which can be 
expressed symbolically as 



±c 

i=l J ° 

(; 

/ - \ T 

far \ 


( 1 
\di>i) 

6\j) 

f dL 

_i’ 

l3ut(4,*)J 

f dL 

J 


5u« + 


dLj 

#u' 


K + 


N 


dxi+^2 


dL 


i T 


r di y 

.<9V>i(4,0. 


<5u»(£i , t) 
(31) 


Before proceeding with the derivation of La- 
grange’s equations by means of the extended Hamil- 
ton’s principle, Eq. (14), it is advisable to identify a 
set of generalised coordinates capable of describing the 
motion of the system fully. From Eqs. (3), we con- 
clude that the motion of only one of the points Oi is 
independent. We choose this point as 0\ t so that we 
retain only R o i(0 for inclusion in the set of general- 
ized coordinates. On the other hand, because 0* repre- 
sent hinge points, the rigid-body rotation vectors 0<(t) 
(i = 1,2 ,...,jV) are all independent. Similarly, the 
nonzero components of the elastic displacement and ro- 
tation vectors, u^(Xg,t) and V\(*i,0 (* = 1,2,.. . , N)> 


where L, = T» - VS is the Lagrangian density for 
body i. Moreover, (dL/d q) T represents the row ma- 
trix [dL/dqi dL/dq 2 “ dL/dqN R ]> etc., where Nr is 
the total number of independent rigid-body degrees of 
freedom. Consistent with the generalized coordinates 
used, the virtual work has the form 

* rti 

SW =Q T iq + V] / (f/" <* u i + m T <*V\) dXi 

i=l J ° 

+ ^[UT^(A,0 + »ffi^(A,t)] (32a) 
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where we write the generalised force vector Q in the 
form 

Q = [F? Mf Mj • ■ • M^] T (32b) 

and note that Fi it a generalised force and Mi . .,M// 
are generalised torques. They can all be related to the 
actuator forces and moments, but we postpone further 
discussion of this subject, and the derivation of specific 
formulas for U< and *4 until later. 

Introducing Eqs. (31) and (32) into Eq. (14), car- 
rying out the usual integrations by parts and recalling 
that the virtual displacements vanish at t = *i, t 3 , we 
have 



dL d dL ’ 

+ £i\t **(*.*) af[au<(A.t). 

+u ‘} T< «‘< 4 '‘) + {a^b) 

-M^] + *‘r < ^' ,) )} i,=o(33) 


Then, invoking the arbitrariness of the virtual displace- 
ments, we obtain the system Lagrange’s equations of 
motion 



where u« and V’, must be such that the equations 



i= 1,2,... ,1V- 1 (35c) 



+ *i 

n=ii 


(1 

' dL ' 

\dt 



dL 


}) T **(4,t) = 0, 

• = 1,2 JV-1 


(35d) 




= 0 


*N=tN 


= 0 


*S=tN 


(35e) 

(35f) 


must be satisfied. Recalling that the body axes 
are embedded in the body at = 0 , we conclude that 
satisfaction of Eqs. (35) is guaranteed if 


Ui(0,*) = 0, *,(0.0 = 0, » — 1, 2, .... JV (36a,b) 


d 

‘ dL 1 

CO 1 
1 

CO 

dt 

.5u<(4,t)j 

d\n(ii,t) 5u<| 



* = 1,2,...,//- : 

d 

‘ dL ' 

dL _ dLi 

dt 


Si /\(4,0 ~ d\l>[ 


+ U <( 
(36c) 


+ •<. 




« = 1,2,...,JV-1 


(36d) 


dl N 

du' N 


= 0, 


dL N 


Wn 


= 0 




(36e,f) 


Equations (34a) represent ordinary differential equa- 
tions for the rigid-body motion and Eqs. (34b) 
and ( 34 c) represent partial differential equations for the 
elastic motions. Moreover, Eqs. (36) are recognized as 
the boundary conditions accompanying the partial dif- 
ferential equations. Although Eqs. (34a), Eqs. (34b), 
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(36a), (36c) and (36e) on the one hand and Eqs. (34c), 
(36b), (36d) and (36f) on the other hand have the ap- 
pearance of independent sets of equations, they are in 
fact simultaneous. They constitute a hybrid (ordinary 
and partial) aet of differential equations governing the 
motion of the multibody system shown in Fig. 1. 


4. Lagrange’s Equations for Flexible Multi- 
body Systems in Terms of Quasi-Coordi- 
nates 

Equations (34) seem very simple, but they axe not. 
The reason for this is that the kinetic energy is only an 
implicit function of q and q and not an explicit one. 
The kinetic energy is an explicit function of V<* and 
, which are commonly known as derivatives of quasi- 
coordinates (Ref. 33). Actually, the kinetic energy is an 
explicit function of n<, but is related directly to u;*, 
as can be seen from Eq. (13). As shown in Ref. 32 for 
a single flexible body, hybrid Lagrange's equations of 
motion in terms of quasi-coordinates are considerably 
simpler than the standard Lagrange's equations. We 
propose to show in this paper that the same is true for 
multibodies. 

Recalling definition (29) of the rigid-body displace- 
ment vector q(t), we can rewrite Eq. (34a) in the more 
detailed form 


dt 



=Fi (37a) 

=M<, » = 1,2, ... (37b) 


The vectors R 0 i, R 0 i and Fi are in terms of compo- 
nents along the inertial axes XYZ. Moreover, the com- 
ponents of the symbolic vector 0* represent rotations 
about nonorthogonal axes leading from to 
and the components of are associated moments. An 
example of such rotations are Euler's angles (Ref. 33). 
As the quasi- velocity counterpart of the generalised ve- 
locity vector q(t), we choose 

w = [Vj, uj uj ■ • • v T N ] T (38) 


and we note that w does not equal the time derivative 
q of the displacements. We also note that every 
three-dimensional vector entering into w is in terms 
of the corresponding orthogonal body axes x^Zi. The 
relation between the velocity vector V<,x in terms of 
body axes and the velocity vector R<,i in terms of 
inertial axes is simply 

V.i = CjR ol (39) 


ia term* of body axes and the Eulerian-type velocities 
8i can be written as 


Vi — D{6i, i — 1, 2 N (40) 

where D{ ia a given tranaform&tion matrix (Ref. 33). 
Equation* (39) and (40) and their reciprocal relations 
can be expressed in the compact form 

w = A T (q)q, q = 5(q)w (41a,b) 

where 

A = block-diag[Cf D* D% • * • d£] (42a) 

B = block-diag[Cf Df 1 Dj 1 * • D^ 1 ] (42b) 

Equations (37) postulate a Lagrangian in terms of 
generalised coordinates and velocities, Eq. (30), when 
in fact the Lagrangian defined by Eqs. (15), (16), (17), 
(25) and (26) is in terms of generalised coordinates and 
quasi- velocities. To distinguish between the two forms, 
we define 


L* =I>*[q,w,Uj,u^,Ut, 

w(4, 0> w(4, 0* V><(4, 0» V’iCA, *)] (43) 


We propose to obtain Lagrange’s equations in 
terms of quasi-coordinates by transforming Eqs. (37). 
To this end, we use the chain rule for derivatives with 
respect to vectors and consider Eq. (39) to obtain 


6L djCi Rot) 7 dL* _ ~ dL' 
dRai ~ dRoi 5V.i ~ Gl 5V 0l 
dL d£_ 

5Roi dRoi 


(44a) 

(44b) 


But, it is shown in the Appendix that the matrix of 
direction cosines C* and quasi-velocity vector u ^ satisfy 
the relation 

Ci = QjCi (45) 

so that differentiating Eq. (44a) with respect to time, 
we have 



Then, inserting Eqs. (44b) and (46) into Eq. (37a) 
and premultiplying by C\ t we obtain the translational 
Lagrange's equations in terms of quasi-coordinates 



where C\ is the matrix of direction cosines first intro- 
duced in Sec. 2, and that between the velocity vector 
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where 


where 


MJ = £>f T M,-, i = 2,3,...,N 


(56) 


FJ = C x F t (48) 

ia the resultant force acting on body 1 in terms of body- 
axes components. 

As far as the rotational motion is concerned, we 
consider first the equations for body 1. Using the chain 
rule for derivatives with respect to vectors once again 
and using Eq. (40), we obtain 


d0i 30 1 du>i 1 3u i 

3L _ 3L* 0(CiR,i) t 3L* 3(Di0i) t 3L* 

30 1 “ 30 i + 30 i av 0l + 30 1 3 ui 

(49b) 

Moreover, Eq. (A-29) from the Appendix, with a re- 
placed by R*i yields the relation 


«(C 1 R.,) r „r,-, 

W ~ ' 


and Eq. (A-27) shows that 

\T d(Di0if 


D[ = 


+ i 


de Y 


(50) 


(51) 


Hence, using Eqs. (49)-(51), we can write 


(9L\ 

3L 

- dT - 

d(Di0if' 

dr 

\30iJ 

30 1 

— u \ 

30 1 


_<r d 

I '31 L* 


3L* 3L * 

+ d ‘5 


) + 1 V ° l 3V 0 1 30 1 

°*li( 

3L r\ 

dui) 

- dir 

+ Vol 3\ 0 i 

. er l 

•+«i - — 
uWj 

<«> 


Inserting Eq. (52) into Eq. (37b) and premultiplying 
the result by , where the superscript -T denotes 
the inverse of the transposed matrix, we obtain the 
rotational Lagrange's equations for the first body in 
terms of quasi-coordinates 



where 

MJ = 0f r M! (54) 


is the resultant torque acting on body 1 in terms of 
body-axes components. The equations of motion for 
the remaining bodies can be obtained in the same man- 
ner, except that V* (i = 2, 3, . . . , N) are not indepen- 
dent, as can be concluded from Eqs. (10). Hence, from 
Eq. (53), the remaining rotational Lagrange's equations 
in terms of quasi-coordinates are 


d (dL'\ m 8L • 0 _ t S2/ 

dt +Wi du>i Di 30i ~ M< ’ 

i = 2,3,...,// 


Equation! (47), (53) and (55) can be cast in a (in- 
gle matrix equation. Indeed, recalling Eqa. (29), (38), 
(41b) and (42b), the rigid-body Lagrange’s equation! 
of motion in term of quasi-coordinates can be written 
in the compact form 



where the asterisk in L* wan dropped for convenience. 
Moreover, 


f wi 

V„i 


0 

0 


0 

0 

<3, 


0 ■ 

0 

0 


L o o o ■■■ Qn J 


and 


(58) 


Q* = fl r Q = [FI T MI T M; T -.M^ r ] T (59) 


The hybrid set of equations of motion is completed by 
adjoining to Eq. (57) the partial differential equations 
for the elastic motions, Eqs. (34b) and (34c), and the 
associated boundary conditions, Eqs. (36). 


5. Explicit Hybrid Equations of Motion for 
Flexible Multibody Systems 

Using Eqs. (16) and (17), we can write the kinetic 
energy in the form 


'=£/ *< d*i 

i N r li 

= - V' [miV£V oi +n£j«n r < + / Piiijiiidxi 

+ f jdifri dx. + 2V^(S^* flr« + / p.u. dx{) 

Jo Jo 

+ 2 n£ J*' (h* + dx<j 


(60) 


and we observe that T does not depend explicitly on 
the quasi-velocities V a , and Ui (i = 1, 2, . . . , N), but 
on V a< and ft ri (i = l,2,...,Af). To resolve this 
inconvenience, we make use of the discrete step function 
fi, defined by 


(55) 


7 


_ JO, iff = -1,-2, -3,... 
7< " 1 1, if » = 0, 1,2,3, ... 


(61) 



and then make repeated use of Eqs. (10) and (13) to 
establish the relations 


N 

Or < + 7i-j-i^,j (*j . 01 ( 62a ) 

i = i 

N 

Vo< =C? 1 V ol + [fiJOri + v,-(/y,0] 

JV 

=Qv 0l + 53 {r y w> + ru + in„(o,0 

;=i 

+7<-i-iC y v J (/ ;i t)} (62b) 

Ori =53 C?i [?i W; + Ti - 1 0 y (^ , t )] + ^ ( 63& ) 
J=i 

N 

V 0i =c*,v oX + 53 {r yW/ + r iJ+1 n y (^,t) 

+7 < -;_iC7 i v ; (/ ;) 0} + dvi (63b) 

N 

6®U =£c? [Tw^i + 7<-j-i^(4,0] (63c) 

;=i 

AT 

mv =c; 1 «r ; 1 + 53 [r y w; + r iJ+ i^,(4,0 

; =1 

+ y i . i .iq i Su i (t i ,t)] (63d) 

in which C* ; - is simply the matrix of direction cosines 
of axes XiyiZi with respect to axes XjyjZj t defined for 
all indices i } j between 1 and N } and consequently 

<*- II Ct % 1 <j<i<N (64*) 
*=;'+! 


C'i = I, 1 < i < N 


(64b) 


(C- i ) T = q <) CT^Cl^Cti, 1 <iJ,k<N 

(64c, d) 

The other quantities appearing explicitly or implicitly 
in Eqs. (62) and (63) are given by 


(65a) 

(65b) 


Uei — Uyt(£j,t) Uji(£i,t)] 

r# =E c ?»«2i c u 

*=; 

N 

dn* = C*j [t »-;W; + 7 »-;-ift«; (lj , <)] (65c) 

i-i 

i- 1 

dv< =c7iV„i + 53 {<%[«£ n ri + vj(<j,0] 

i=i 

+CJ [ft, i »,(<„l) + S5<i<v]} («M) 

CJ = (-0,, + CJlC,) CJ («5e) 


We also note that depends only on 0k, for 
min(i,.j) < * < max(i,i), and on ^ k (*k,t). for values 
of k satisfying min(t, j) < k < max(»,j). Hence, using 
Eqs. (A-29) and (A- 30), we can derive the relations 


d -^f- = (Ji-U - yi-kWlci; sq { (66a) 


provided a does not depend on 0*, and 

- 7<-*-i (66b) 

provided a does not depend on ^ k (l*>*)' Some other 
relations that will prove useful are as follows: 


dnl 


= 0 


SR 0 i 

7=* = D l - yi-k)Ckjbi-j*i 


(67a) 


j=i 


dnl 


= o 


+ 7 < _,_ 1 n. i (4,t)]q < (67b) 

w (67c) 

auk(/*,t) v 

onT N 

+ y i - j - l n. i (t i ,t))q i (67d) 

^=o,g= 7i .,cr. ( 67«,f) 

dn ^~ = 0, A = (67g,h) 


dvk(/s,t) ’ wi*(4,t) 


5V T 

= 0 

3R«i 

^- = D k (7i-k-7i-k)C;,V 0l C' u 


(68a) 


|j|(7 i -k-7i-0C;y [^ n -7+ v ;(< 


as— £ k» 


0u*(4k,O 

0V£ 


(68b) 

(68c) 


;=i 


+ 


= £*(**,<) 


— 7i— *— 


^{(^-k-l -7i-k-i)<% fanrj+VjiCt) 
J=1 l 


8 


dvl 


dVol 


+ Ekl(yti,t ^di/} k (t k 

dvl 

du k ~ a 




0 % _ r*+ v w oi — r»T 
— ^1»» - Aii 

d vjj _ _ d\l 


dv k (lk,t) 


— yi-u-xCl, 


50 ,*(/*, f) 


(68d) 

(68e,f) 
= i?* +1 (68g,h) 


Then, using the chain rule for vectors when needed, we 
obtain the momenta 


TV 


m dL L 

ffVoi wl '5Voj 


_ dL 

PVol =^7— - Z^ c 
1=1 


d - QL - /"r T — it c* ^ 

Ph " ~a W; ~ ^ V ,J av oi + 7, - ;Gi > an ri } 


(69a) 

(69b) 


where 

dLi 

av„ 

5Li 


= m*V 0 j + Sf flri 4” / Pi d* dx^ 
n JO 


(70a) 


=Jt<n r< + S<V a< + ^ ’ (!<!»< + J^) dxi (70b) 
For future reference, we also indicate that 

Tt{w~) + 

+ / piUidxi + d t vi (71a) 

Jo 

i(^“) = J "«rt + $|V oi 

+ j (Ijiii + dzj + d,ni (71b) 


where 


d|Vi — f PiUidXi (72a) 

*/o 

dtfk = jti^ri ““ K>» f Pid%dZ{ (72k) 

Jo 

jfi =j Pi{M*iZi + +(*& + UifiT] dZi (72c) 

and 

Pv«i = "*») + E 

+r i -/C7 j sTc;)] *, 

tv r tv 

Erw-imi^l v ; (^-,t) 


4=1 


+E 

i=i 


(73a) 


TV 


o 1 


tv r n 

+ y: 

; = 1 L*=l 

+T*-y - 1 C*i C. r / ) ] 0,y(^, t) 

+ E J Q Pi^ij 

+ ^[^u^- + C t u (rn i d Vi 

+5^ dm + dtv»J j 

P*; = I £ (^r£ + Td-iCjSi) C?] 

Li=l 

N f N 

+ EiE[( m<r 5 + 7wC7 i 5 < )r < * 
k*l l »=1 

+ 7«-* (rJjST + 7i-}CijJti) CJj] |<*»* 

n r n 

*=i Li=i 

Cfc] v*(/*,0 

TV r TV 

+ £ I E [ ("* r 5 + *->c&$) r<,*+i 

h=l l i=l 

*»-*-! (rZST+n-jCSJu) c?*] 


+ 


TV 


yii . .. 

n r 

/ Jei^idZi 

+E 

Jo 

»=i L 


+ Ji-jCij q^~ + rfjdtvi + TwC^dtfu + ( m »r?y 

+Ti— yCJy5») dv» + ^ 7i— j‘Cy •Tt*) do» 

(73b) 

We also define equivalent forcea and momenta 

f;i= c i^ = ° ( 74a ) 


L aR*, 


m;, 


-n~ T QIl — n~ T V' f ai> dL 

~ j d6j - > ^ \ as. av.. aa. an.. 


^ V aa ; - av oi + aa ; an 

(74b) 

and the remaining pertinent terma 


dL A avr dL t 

tohViS) ^auy^.t) av« 


(75a) 
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dy^t) kdvMMV* 1 > 

ql _^( 0v£ «L an£ 5L \ 
tyfc.O “fij Wi«ht)«Vw + ty(4,<) «<W 

(75c) 

aL avt ai an?; sl \ 

«M4.0 “£? Van.,^, 0 av #< an. ; (^,t) an*) 

(75d) 

in which some of the partial derivatives are given by 
Eqa. (87). 

Finally, adjoining the kinematic relations ex- 
pressed by Eqs. (9), (11), (39) and (40) and inserting 
Eqs. (68)-(70) into Eqs. (34b), (34c) and (57), we obtain 
the hybrid state equations in terms of quasi-coordinates 

R,i=C?V ol , 9i=Dr l Ui , 1*1,2,...,* 

(76a, b) 

Ui(xj,t) = ▼<(*<, t), ^j(*i,t) = n e< (xi,t), 

i = l,2,...,* (76c, d) 

PVoi = — wipvoi + F* (76e) 

P «1 = -v.iPvoi - W 1 P «1 + m;, + (76f) 

Pwi = -WiPw* + M*< + MJ<, i = 2, 3, .... * (76g) 

"b 4" — Ujilirai — 2f2 rt ,l', l + Q r j« V« t i 

firs, Vojt + (f^r.i fl rl i)t*yi 

+ flry»fV.M u *»] — [kyiGii4,(u y i — V’*,)] — fyi (76b) 


dLj 

au;. 


. . **(4,0/ 

***** 


*= 1,2, ...,* — 1 


(77c) 


BLi 

d+\ 


*i=r* 


_ n r dL i _ \ _ * 

\dt |an.,<(4,oJ **(*,0/ 

» — 1,2, ...,*— 1 (77d) 


dL N 


du' N 


= 0, 


dLs 






= 0 (77e,f) 


*W=^N 


and the generalised forces and torques are given by 
Ft=Ecr i F; i (78a) 

4=1 

m* =X)(rf 1 F; i + cj i M; < ) (78b) 


i=l 




M,’ =MV + £ ( r f‘ F *i + 7,-iC-^MV) , 

;=i 

» = 2,3,...,JV (78c) 

Ui=£r;.4-iC7 i F;., » = (78d) 

;=i 

*i =£ (ifr+iF;, + , 

; = 1 

»= 1,2,...,//-! (78e) 


where we have made use of Eqs. (27), (32a) and (63c, d). 


Pi[v M % + — Zgflryft + + 2flrri v yi + ^rzi^oyi 

— Q r ygV»t + *• “ (^r#» 4* ^ryi) 14 ** 

4* flryiflrji^vi] “ + ^y»)] = f»i (76i) 

J ziaii&ami 4" (irai) “ (^«*Gi/**^ c ,) = (76j) 

Jyiyi(Q*yi 4“ flryi) 4* ^aidiAi(u g ^ + ^yi) — (Ej/y*t/iyj) 

= nty(76k) 


4* “ 4y»G|Aj(u^ — l/>jj) 

= m J *(76k) 

The associated boundary conditions, Eqs. (36), are 
given by 

114 ( 0,0 = 0, *,(<),*) = 0, 1 = 1,2,...,// (77a,b) 


6. Summary and Conclusions 

In recent years, there has been an increasing interest in 
deriving the equations of motion for flexible multibody 
systems by treating the mass and stiffness of the bod* 
ies as distributed parameters. The equations of mo- 
tion are generally derived by means of the extended 
Hamilton’s principle, leading to a hybrid set of equa- 
tions, where hybrid is to be taken in the sense that 
the rigid-body translations and rotations of the bod- 
ies are described by ordinary differential equations and 
the elastic motions are described by partial differential 
equations with appropriate boundary conditions. In 
earlier investigations, the rigid-body rotations were de- 
scribed by Eulerian-type angles, which tend to compli- 
cate unduly the equations of motion, unless the motion 
remains planar. 

This paper presents a mathematical formulation 
for flexible multibodies in terms of quasi- coordinates, 
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which permit* the derivation of the equations for gen- 
eral rigid-body motion* with consider ably more ease 
than be using Eulerian-type angles. As an added fea- 
ture, the equations for the elastic motions include rota- 
tory inertia and thegr damnation effects. The equa- 
tions of motion are- npt in state form, making them 
suitable for control 
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Appendix 
i. Derivative rules 


If A = [Aij] is an m x n matrix, then we define the 
partial derivative of A with respect to a scalar r to be 
the m x n matrix dA/dr as [0Ay/0r]. If A is a function 
of time t, then the derivative of A with respect to t is 
denoted by A = dA/dt = [dA^y/dt]. Let B = [fly] be 
nnMxN matrix. Then, the derivative of a matrix with 
respect to a matrix, 0A/6B, is the mM x nN matrix 
defined by 


9A 

dA 

dA -j 

dB 11 

dBu 

dB in 

dA 

dA 

dA 

dBn 

dBu 

8 Bin 

dA 

dA 

dA 

'■dBu 1 

dB U 2 

dBuN 


(A — 1) 


Furthermore, let L be a scalar and f = [f\ ■ • • f m ] T , 
q = [(j... y n ] T , s = [*!••• Sr] 5 * be column matrices. 


Then dL/dq is a column matrix, df r /dq is an n x m 
matrix and df/dq r = (df r /dq) T . The chain rules for 
differentiation have the form 

dP_ d<? dF 9t _ dt dq 

8% 3s 3q ° r du T ~ 3q T 8u T 

8L _ 9q T dL dL dL <9q 
dn ~ d* 0q ° F d» T ~ dq T 8u T 
Moreover, 

A dr et . 
t= w* 

%rT*- = A o, = 

dq T oq 

dq 

provided A does not depend on q. 
iL Proper orthogonal matrices 

Throughout this paper, we encounter proper orthogo- 
nal matrices C, which are functions of three indepen- 
dent coordinates 9 = [0 X 0 a PsF- These matrices can 
be identified as matrices of direction cosines of one co- 
ordinate system with corresponding unit vectors 

bi, b 3 , bs, with respect to another coordinate system 
*1*2*3, with corresponding unit vectors ni, 0|, ns. 
Hence, letting C = [Cy], the entries Cy can ie ex- 
pressed as 


(A -2) 
(A -3) 

(A -4) 
(A -5) 
(A -6) 


Cij = bj • ny , i,j — 1,2,3 
which implies that 


(A -10) 


a 3 

ny = y*l(ny • b*)b* = Clyb*. j = 1, 2, 3 

At this point we wish to establish a relation between 
the body axes components of the angular velocity w 
of coordinate system (1(2(3 with respect to coordinate 
system *1*3x3 and the time derivative of Cy with 
respect to coordinate system *1*3*3- First, recall 
(Ref. 33) that u is uniquely characterised by 


K, = u x b<. i = 1.2,3 


(A - 121 


where in this case the "dot” requires holding ni , n a , ns 
constant. Then, taking the time derivative of Eq. (A- 
10), using Eqs. (A-ll) and (A-12), and some identity 
involving scalar and vector products, we obtain 


Cy =b* • ny = (<•> x b<) • ny = (b< x ny) • u 
3 3 

=(b 4 x y^Ctybt) • v = ^2Ckj(bi x b*) • u 

(A - 13) 
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kml 
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Now we observe that (b* x b») • w, where i, k s 1, 2, 3, 
ere merely the entries of the 3 x 3 matrix 



where [wi wj wa] T are the fxfjfs components of w, 
and we have used the fact that bx, b?, bj form a right- 
handed set of unit vectors. Inserting Eqs. (A- 14) into 
Eq. (A-13), we obtain 

3 

Cij = £<**<?*, (A -15) 

tel 

which can be expressed in the matrix form 

C = Z T C (A - 16) 


The relationship between w and 8 has the form 


u = D{9)9 


(A -17) 


We now propose to derive tome relation between D 
and C. In the first place, taking the partial derivative 
of CC 7 * = I with respect to we obtain 

c°$L + - c^- + (c^\ - 0 

c Mi + Mi ~ C Mi +V Mi ) °* 

i= 1,2,3 (A -18) 

from which we conclude that the 3x3 matrix 
C(dC T /Mi) is skew symmetric. We denote the ma- 
trix by 


§i = c W ' i = 1 ' 2 ’ 3 


(A -19) 


where Si is obtained from the column matrix S< = 
[5x< S% in the usual manner. We now calculate 
the time derivative of C? in the farm 

tel * tel V * ' tel 

=C? S<*<) = ((ft S3 Sa] 9) = C?(S9) 

(A -20) 

Comparing Eqs. (A- 16), (A- 17) and (A-19), we con- 
clude that 


5= [Si S, S 8 ] = D 


(A -21) 


Equation (A-20) relates C and D in an implicit manner. 

Next, we wish to derive an ciprasskwi far D. 
Ttking the partial derivative of Eq. (A- 18) with respect 
to 8j and replacing Si by Di, we obtain 


-(«©(«© 


_ PC? nr* . „d*C? * * „PC? 

+ ° 99 j Mi ~ + ° 99 )Mi ~ + C M } Mi 

(A -22) 

Interchanging i and j in Eq. (A-22), we have 


9D, - - PC? 

99i m DiDj + °99 i M j 


(A -23) 


(A -25) 


Then, subtracting Eq. (A-23) from Eq. (A-22), we can 
write 

W"^ =A6 '" 6 ' a=( ®‘ 5 ' ) ( ' 4 ' 24) 

which implies that 

This formula can be used in turn to derive an expression 
for D. First, we recall Eq. (A-17) and write 

= t; ($' ,+6,D, * , i 

= ^ +6 ‘"“^ i+ "’' D ‘ (A - M) 
This implies that 

+ (A-2 7) 

DJ 99 

Next, we consider the partial derivative of (Ca) T 
with respect to 9, where a does not depend on 9. First, 
we recall Eqs. (A-19) and (A-21) and write 

^=* r w=* Ttfr ( c w)-( c -) r °- 

= - (Cm) T £>? = (CaD<) T = -(D,Ca) r 

=Df (Ca) T (A - 28) 

which implies that 

d(Cm) T -Df(^) 

— Dj (Ca) = -D T (Ca) (A -29) 

L-Dj(Ca). 

The companion formula 

= (A — 30) 

vW 

can be derived in a similar manner. 
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